COVID Statistical Analysis
Hongyin Lai, Yuan Li, Tara Prezioso, Alison Rector, Jeffrey Brennan, Swaminathan Kumar, Jose-Miguel Yamal
set.seed(1)# Data manipulation & cleaning
library(tidyverse)
library(reshape2)
library(nlme) # gapply
library(lubridate) # date cleaning
# Viz
library(gridExtra)
library(ggmap)
library(maps)
library(mapdata)
# Modeling & stats
library(mgcv) # gam
library(R0)
library(Kendall) # Mann-Kendall
# Time series & forecasting
library(forecast)
library(zoo)
library(astsa)
library(fpp2)TODO
- combine time series predictions into one generalized function & improve filenames/paths
RT ANALYSIS
#read in full DSHS data with county and TSA information
# TOOD: rename
full.dshs <- read.csv("combined-datasets/county.csv")
tsa.dshs <- read.csv("combined-datasets/tsa.csv")
phr.df <- read.csv("combined-datasets/phr.csv")rt.data.clean<-function(covid.data) {
# select relevant metadata & DSHS population estimates
rt.cols <- c("County", "Date", "TSA_Name", "TSA", "PHR", "PHR_Name",
"Metro_Area", "Cases_Daily", "Population_DSHS")
covid.data <- covid.data[, names(covid.data) %in% rt.cols]
# set correct types
covid.data$Date <- as.Date(covid.data$Date)
covid.data$Cases_Daily <- as.numeric(covid.data$Cases_Daily)
# TODO: decide on how to handle drops in cumulative counts -> negative daily cases
covid.data <- covid.data %>%
mutate(Cases_Daily = ifelse(Cases_Daily < 0, 0, Cases_Daily))
return(covid.data)
}#separate county, metro, TSA and state data into separate dataframes
rt.data.organize <- function(mycounty, mytsa, myphr){
### COUNTY ###
cleaned.county <- rt.data.clean(mycounty)
county.df <- cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area, -TSA)
### TSA ###
TSA.df <- rt.data.clean(mytsa)
### PHR ###
PHR.df <- rt.data.clean(myphr)
### METRO ###
# calculate daily cases and population at metro level, drop repeated rows
metro.temp <- cleaned.county %>%
group_by(Metro_Area,Date) %>%
mutate(metro_cases_daily=sum(Cases_Daily, na.rm=TRUE)) %>%
mutate(metro_pop_DSHS=sum(Population_DSHS, na.rm=TRUE)) %>%
dplyr::select(Date, Metro_Area, metro_cases_daily, metro_pop_DSHS) %>%
distinct()
metro.df <- data.frame(Date = metro.temp$Date,
Metro_Area = metro.temp$Metro_Area,
Cases_Daily = metro.temp$metro_cases_daily,
Population_DSHS = metro.temp$metro_pop_DSHS)
### STATE ###
# calculate daily cases and population at state level, drop repeated rows
state.temp <- TSA.df %>%
group_by(Date) %>%
mutate(state_cases_daily=sum(Cases_Daily, na.rm=TRUE)) %>%
mutate(state_pop_DSHS=sum(Population_DSHS, na.rm=TRUE)) %>%
dplyr::select(Date, state_cases_daily, state_pop_DSHS) %>%
distinct()
state.df <- data.frame(Date = state.temp$Date,
Cases_Daily = state.temp$state_cases_daily,
Population_DSHS = state.temp$state_pop_DSHS)
### OUTPUT ###
rt.df.out <- list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df, phr = PHR.df)
return(rt.df.out)
}rt.df.extraction <- function(Rt.estimate.output) {
# extract r0 estimate values into dataframe
rt.df <- setNames(stack(Rt.estimate.output$estimates$TD$R)[2:1], c('Date', 'Rt'))
rt.df$Date <- as.Date(rt.df$Date)
# get 95% CI
CI.lower.list <- Rt.estimate.output$estimates$TD$conf.int$lower
CI.upper.list <- Rt.estimate.output$estimates$TD$conf.int$upper
#use unlist function to format as vector
CI.lower <- unlist(CI.lower.list, recursive = TRUE, use.names = TRUE)
CI.upper <- unlist(CI.upper.list, recursive = TRUE, use.names = TRUE)
rt.df$lower <- CI.lower
rt.df$upper <- CI.upper
rt.df <- rt.df %>%
mutate(lower = replace(lower, Rt == 0, NA)) %>%
mutate(upper = replace(upper, Rt == 0, NA)) %>%
mutate(Rt = replace(Rt, Rt == 0, NA))
return(rt.df)
}#Function to generate Rt estimates, need to read in data frame and population size
covid.rt <- function(mydata) {
### DECLARE VALS ###
#set generation time
#Tapiwa, Ganyani "Esimating the gen interval for Covid-19":
# LOGNORMAL OPS
# gen.time<-generation.time("lognormal", c(4.0, 2.9))
# gen.time<-generation.time("lognormal", c(4.7,2.9)) #Nishiura
# GAMMA OPS
# gen.time<-generation.time("gamma", c(5.2, 1.72)) #Singapore
# gen.time<-generation.time("gamma", c(3.95, 1.51)) #Tianjin
gen.time <- generation.time("gamma", c(3.96, 4.75))
print(as.character(mydata[1,2]))
# TODO: consider removing this and handling na cases directly
#change na values to 0
mydata <- mydata %>% replace(is.na(.),0)
sum.daily.cases <- sum(mydata$Cases_Daily)
pop.DSHS <- mydata$Population_DSHS[1]
#Get 7 day moving average of daily cases
mydata$MA_7day<-rollmean(mydata$Cases_Daily, k=7, na.pad=TRUE, align='right')
#change na values to 0
mydata<-mydata%>%replace(is.na(.),0)
#create a vector of new cases 7 day moving average
mydata.new<-pull(mydata, MA_7day)
#mydata.new <- pull(mydata, Cases_Daily)
# get dates as vectors
date.vector <- pull(mydata, Date)
#create a named numerical vector using the date.vector as names of new cases
#Note: this is needed to run R0 package function estimate.R()
names(mydata.new) <- c(date.vector)
#Find min.date, max.date and row number of max.date
min.date <- min(mydata$Date)
max.date <- max(mydata$Date)
#get row number of March 15 and first nonzero entry
#NOTE: for 7 day moving average start March 15, for daily start March 9
#find max row between the two (this will be beginning of rt data used)
march15.row <- which(mydata$Date=="2020-03-15")
first.nonzero <- min(which(mydata$Cases_Daily>0))
last.nonzero <- max(which(mydata$Cases_Daily>0))
minrow <- max(march15.row, first.nonzero)
### R0 ESTIMATION ###
# TODO: establish better criteria for # of required daily cases
# mean daily cases; average last x number of days; % of region pop etc.
if(!is.na(minrow) & sum.daily.cases > 100) {
# declare limits for rt estimation (first nonzero date & last nonzero date)
i <- minrow
j <- as.integer(min(last.nonzero, max.date))
#reduce the vector to be rows from min date (March 9 or first nonzero case) to current date
mydata.newest <- mydata.new[i:j]
rt.DSHS <- estimate.R(mydata.newest,
gen.time,
begin=as.integer(1),
end=length(mydata.newest),
methods=c("TD"),
pop.size=pop.DSHS,
nsim=1000)
rt.DSHS.df <- rt.df.extraction(rt.DSHS)
} else {
# error catch small regions & na values for minrow
rt.DSHS.df <- data.frame(Date = as.Date(mydata$Date),
Rt = rep(NA, length(mydata$Date)),
lower = rep(NA, length(mydata$Date)),
upper = rep(NA, length(mydata$Date)))
}
return(rt.DSHS.df)
}County
#apply the covid.rt function to the county.datily data by county
county.rt.output <- nlme::gapply(county.daily, FUN = covid.rt, groups = county.daily$County)## [1] "Anderson"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Andrews"## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.142857142857143, `2020-04-05` =
## 0.857142857142857, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.142857142857143, `2020-04-05` =
## 0.857142857142857, : Using initial incidence as initial number of cases.## [1] "Angelina"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Aransas"## Warning in est.R0.TD(epid = c(`2020-04-05` = 0.142857142857143, `2020-04-06` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-05` = 0.142857142857143, `2020-04-06` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Archer"
## [1] "Armstrong"
## [1] "Atascosa"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Austin"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Bailey"## Warning in est.R0.TD(epid = c(`2020-05-05` = 0.142857142857143, `2020-05-06` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-05-05` = 0.142857142857143, `2020-05-06` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Bandera"
## [1] "Bastrop"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Baylor"
## [1] "Bee"## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Bell"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Bexar"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Blanco"
## [1] "Borden"## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf## [1] "Bosque"
## [1] "Bowie"## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Brazoria"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Brazos"## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Brewster"## Warning in est.R0.TD(epid = c(`2020-05-01` = 0.142857142857143, `2020-05-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-05-01` = 0.142857142857143, `2020-05-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Briscoe"
## [1] "Brooks"
## [1] "Brown"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Burleson"## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Burnet"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Caldwell"## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Calhoun"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Callahan"
## [1] "Cameron"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Camp"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Carson"
## [1] "Cass"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Castro"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Chambers"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Cherokee"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Childress"
## [1] "Clay"
## [1] "Cochran"
## [1] "Coke"
## [1] "Coleman"
## [1] "Collin"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Using initial incidence as initial number of cases.## [1] "Collingsworth"
## [1] "Colorado"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.285714285714286, `2020-04-02` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.285714285714286, `2020-04-02` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Comal"## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.428571428571429, `2020-03-23` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.428571428571429, `2020-03-23` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Comanche"
## [1] "Concho"
## [1] "Cooke"## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Coryell"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Cottle"
## [1] "Crane"
## [1] "Crockett"## Warning in est.R0.TD(epid = c(`2020-06-05` = 0.142857142857143, `2020-06-06` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-06-05` = 0.142857142857143, `2020-06-06` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Crosby"
## [1] "Culberson"
## [1] "Dallam"## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Dallas"## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 1.14285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 1.14285714285714, : Using initial incidence as initial number of cases.## [1] "Dawson"## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Deaf Smith"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Delta"
## [1] "Denton"## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "DeWitt"## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.142857142857143, `2020-03-20` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.142857142857143, `2020-03-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Dickens"
## [1] "Dimmit"
## [1] "Donley"
## [1] "Duval"## Warning in est.R0.TD(epid = c(`2020-04-14` = 0.142857142857143, `2020-04-15` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-14` = 0.142857142857143, `2020-04-15` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Eastland"
## [1] "Ector"## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Edwards"
## [1] "El Paso"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Ellis"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Erath"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Falls"
## [1] "Fannin"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Fayette"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Fisher"
## [1] "Floyd"
## [1] "Foard"
## [1] "Fort Bend"## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Franklin"
## [1] "Freestone"## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Frio"## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Gaines"
## [1] "Galveston"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Garza"
## [1] "Gillespie"## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Glasscock"
## [1] "Goliad"
## [1] "Gonzales"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Gray"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Grayson"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Gregg"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Grimes"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Guadalupe"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.571428571428571, `2020-03-26` =
## 1.14285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.571428571428571, `2020-03-26` =
## 1.14285714285714, : Using initial incidence as initial number of cases.## [1] "Hale"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hall"
## [1] "Hamilton"
## [1] "Hansford"
## [1] "Hardeman"
## [1] "Hardin"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Harris"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Using initial incidence as initial number of cases.## [1] "Harrison"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hartley"
## [1] "Haskell"
## [1] "Hays"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hemphill"
## [1] "Henderson"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hidalgo"## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.285714285714286, `2020-03-24` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.285714285714286, `2020-03-24` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Hill"## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hockley"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Hood"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Hopkins"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Houston"## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Howard"## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hudspeth"
## [1] "Hunt"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Hutchinson"## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.285714285714286, `2020-04-05` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.285714285714286, `2020-04-05` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Irion"
## [1] "Jack"
## [1] "Jackson"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Jasper"
## [1] "Jeff Davis"
## [1] "Jefferson"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 1, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 1, : Using initial incidence as initial number of cases.## [1] "Jim Hogg"
## [1] "Jim Wells"## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Johnson"## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Jones"## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Karnes"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Kaufman"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Kendall"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Kenedy"
## [1] "Kent"
## [1] "Kerr"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Kimble"
## [1] "King"## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf## [1] "Kinney"
## [1] "Kleberg"## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Knox"
## [1] "La Salle"## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Lamar"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Lamb"## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Lampasas"
## [1] "Lavaca"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Lee"## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Leon"## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Liberty"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Limestone"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Lipscomb"
## [1] "Live Oak"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Llano"
## [1] "Loving"## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf## [1] "Lubbock"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Lynn"
## [1] "Madison"## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Marion"## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.428571428571429, `2020-04-17` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.428571428571429, `2020-04-17` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Martin"
## [1] "Mason"
## [1] "Matagorda"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Maverick"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "McCulloch"
## [1] "McLennan"## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Using initial incidence as initial number of cases.## [1] "McMullen"
## [1] "Medina"## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Menard"
## [1] "Midland"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Milam"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Mills"
## [1] "Mitchell"
## [1] "Montague"
## [1] "Montgomery"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Moore"## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Morris"
## [1] "Motley"
## [1] "Nacogdoches"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Navarro"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Newton"
## [1] "Nolan"## Warning in est.R0.TD(epid = c(`2020-04-24` = 0.142857142857143, `2020-04-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-24` = 0.142857142857143, `2020-04-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Nueces"## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Ochiltree"
## [1] "Oldham"
## [1] "Orange"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Palo Pinto"## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Panola"## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.428571428571429, `2020-04-03` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.428571428571429, `2020-04-03` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "Parker"## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.142857142857143, `2020-03-24` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.142857142857143, `2020-03-24` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Parmer"## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Pecos"## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Polk"## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Potter"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.285714285714286, `2020-03-22` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.285714285714286, `2020-03-22` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Presidio"
## [1] "Rains"
## [1] "Randall"## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.428571428571429, `2020-03-29` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.428571428571429, `2020-03-29` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Reagan"
## [1] "Real"
## [1] "Red River"## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Reeves"## Warning in est.R0.TD(epid = c(`2020-05-08` = 0.142857142857143, `2020-05-09` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-05-08` = 0.142857142857143, `2020-05-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Refugio"## Warning in est.R0.TD(epid = c(`2020-05-04` = 0.142857142857143, `2020-05-05` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-05-04` = 0.142857142857143, `2020-05-05` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Roberts"
## [1] "Robertson"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Rockwall"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Runnels"
## [1] "Rusk"## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Sabine"
## [1] "San Augustine"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "San Jacinto"## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "San Patricio"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "San Saba"
## [1] "Schleicher"
## [1] "Scurry"## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Shackelford"
## [1] "Shelby"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Sherman"
## [1] "Smith"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Somervell"
## [1] "Starr"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Stephens"
## [1] "Sterling"
## [1] "Stonewall"
## [1] "Sutton"
## [1] "Swisher"
## [1] "Tarrant"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Taylor"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Terrell"
## [1] "Terry"
## [1] "Throckmorton"
## [1] "Titus"## Warning in est.R0.TD(epid = c(`2020-04-03` = 0.142857142857143, `2020-04-04` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-03` = 0.142857142857143, `2020-04-04` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Tom Green"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Travis"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Trinity"## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.285714285714286, `2020-04-05` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.285714285714286, `2020-04-05` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Tyler"
## [1] "Upshur"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Upton"
## [1] "Uvalde"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Val Verde"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "Van Zandt"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Victoria"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.428571428571429, `2020-03-27` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.428571428571429, `2020-03-27` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Walker"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Waller"## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.428571428571429, `2020-03-30` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.428571428571429, `2020-03-30` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Ward"
## [1] "Washington"## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.714285714285714, `2020-03-29` =
## 0.714285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.714285714285714, `2020-03-29` =
## 0.714285714285714, : Using initial incidence as initial number of cases.## [1] "Webb"## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Wharton"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Wheeler"
## [1] "Wichita"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Wilbarger"
## [1] "Willacy"## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Williamson"## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 1.14285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 1.14285714285714, : Using initial incidence as initial number of cases.## [1] "Wilson"## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Winkler"
## [1] "Wise"## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Wood"## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Yoakum"
## [1] "Young"
## [1] "Zapata"## Warning in est.R0.TD(epid = c(`2020-04-07` = 0.142857142857143, `2020-04-08` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-07` = 0.142857142857143, `2020-04-08` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "Zavala"## Warning in est.R0.TD(epid = c(`2020-04-17` = 0.142857142857143, `2020-04-18` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-04-17` = 0.142857142857143, `2020-04-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.#convert list of dataframes to single dataframe with county as id col
county.rt.df <- data.table::rbindlist(county.rt.output, idcol = TRUE)
colnames(county.rt.df)[1] = 'County'Metro
#apply the covid.rt function to the metro.daily data by metro region
metro.rt.output <- nlme::gapply(metro.daily, FUN=covid.rt, groups=metro.daily$Metro_Area)## [1] "Metro"## Warning in est.R0.TD(epid = c(`2020-03-15` = 7, `2020-03-16` = 6, `2020-03-17` =
## 6.57142857142857, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 7, `2020-03-16` = 6, `2020-03-17` =
## 6.57142857142857, : Using initial incidence as initial number of cases.## [1] "Non-Metro"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.#convert list of data frames to single data frame with metro as id col
metro.rt.df <- data.table::rbindlist(metro.rt.output, idcol=TRUE)
colnames(metro.rt.df)[1] = 'Metro_Area'TSA
# plot with shaded confidence intervals
rt.plot<-function(rt.data, caption){
library(ggplot2)
caption<-rt.data$TSA[1]
plot<-ggplot(rt.data, aes(Date, Rt))+
geom_ribbon(aes(ymin=lower, ymax=upper), fill="light grey")+
geom_point(color="blue", size=0.2)+
labs(title=caption)
plot
}#apply the covid.rt function to the TSA.daily data by TSA region
TSA.rt.output <- nlme::gapply(TSA.daily, FUN=covid.rt, groups=TSA.daily$TSA)## [1] "A"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.571428571428571, `2020-03-22` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.571428571428571, `2020-03-22` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "B"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.714285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.714285714285714, : Using initial incidence as initial number of cases.## [1] "C"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "D"## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "E"## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Using initial incidence as initial number of cases.## [1] "F"## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "G"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "H"## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "I"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "J"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.## [1] "K"## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "L"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "M"## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Using initial incidence as initial number of cases.## [1] "N"## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "O"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "P"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "Q"## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.57142857142857, `2020-03-16` =
## 1.71428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.57142857142857, `2020-03-16` =
## 1.71428571428571, : Using initial incidence as initial number of cases.## [1] "R"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.## [1] "S"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "T"## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "U"## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "V"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.#convert list of data frames to single data frame with TSA as variable
TSA.rt.df <- data.table::rbindlist(TSA.rt.output, idcol=TRUE)
# add TSA_Name & fix dates
colnames(TSA.rt.df)[1] <- 'TSA'
TSA.rt.df$Date <- as.Date(TSA.rt.df$Date)
TSA.rt.df <- merge(TSA.rt.df, TSA.daily[, c('TSA', 'TSA_Name', 'Date')], by = c('TSA', 'Date'))
# combine TSA
TSA.rt.df$TSA <- paste0(TSA.rt.df$TSA, ' - ', TSA.rt.df$TSA_Name)
TSA.rt.df$TSA_Name <- NULL
#plot Rt by TSA
nlme::gapply(TSA.rt.df, FUN=rt.plot, groups=TSA.rt.df$TSA)## $`A - Amarillo`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`B - Lubbock`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`C - Wichita Falls`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`D - Abilene`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`E - Dallas/Ft. Worth`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`F - Paris`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`G - Longview/Tyler`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`H - Lufkin`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`I - El Paso`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`J - Midland/Odessa`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`K - San Angelo`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`L - Belton/Killeen`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`M - Waco`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`N - Bryan/College Station`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`O - Austin`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`P - San Antonio`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`Q - Houston`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`R - Galveston`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`S - Victoria`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`T - Laredo`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`U - Corpus Christi`## Warning: Removed 1 rows containing missing values (geom_point).
##
## $`V - Lower Rio Grande Valley`## Warning: Removed 1 rows containing missing values (geom_point).
PHR
phr.rt.output <- nlme::gapply(phr.daily, FUN=covid.rt, groups=phr.daily$PHR)## [1] "1"## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 1.28571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 1.28571428571429, : Using initial incidence as initial number of cases.## [1] "11"## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.## [1] "2/3"## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Using initial incidence as initial number of cases.## [1] "4/5N"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "6/5S"## Warning in est.R0.TD(epid = c(`2020-03-15` = 3, `2020-03-16` =
## 2.14285714285714, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 3, `2020-03-16` =
## 2.14285714285714, : Using initial incidence as initial number of cases.## [1] "7"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "8"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.## [1] "9/10"## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.#convert list of data frames to single data frame with phr as id col
phr.rt.df <- data.table::rbindlist(phr.rt.output, idcol=TRUE)
# add PHR_Name & fix dates
colnames(phr.rt.df)[1] <- 'PHR'
phr.rt.df$Date <- as.Date(phr.rt.df$Date)
phr.rt.df <- merge(phr.rt.df, phr.daily[, c('PHR', 'PHR_Name', 'Date')], by = c('PHR', 'Date'))
# combine phr
phr.rt.df$PHR <- paste0(phr.rt.df$PHR, ' - ', phr.rt.df$PHR_Name)
phr.rt.df$PHR_Name <- NULLState
state.rt.df <- covid.rt(state.daily)## [1] "0"## Warning in est.R0.TD(epid = c(`2020-03-15` = 7.28571428571429, `2020-03-16` =
## 6.28571428571429, : Simulations may take several minutes.## Warning in est.R0.TD(epid = c(`2020-03-15` = 7.28571428571429, `2020-03-16` =
## 6.28571428571429, : Using initial incidence as initial number of cases.export
write.csv(county.rt.df,"statistical-output/rt/county_rt.csv", row.names = FALSE)
write.csv(metro.rt.df, "statistical-output/rt/metro_rt.csv", row.names = FALSE)
write.csv(TSA.rt.df, "statistical-output/rt/tsa_rt.csv", row.names = FALSE)
write.csv(phr.rt.df, "statistical-output/rt/phr_rt.csv", row.names = FALSE)
write.csv(state.rt.df, "statistical-output/rt/state_rt.csv", row.names = FALSE)grouping
colnames(county.rt.df)[1] <- 'Level'
colnames(metro.rt.df)[1] <- 'Level'
colnames(TSA.rt.df)[1] <- 'Level'
colnames(phr.rt.df)[1] <- 'Level'
state.rt.df$Level <- 'Texas'
county.rt.df$Level_Type = 'County'
metro.rt.df$Level_Type = 'Metro'
TSA.rt.df$Level_Type = 'TSA'
phr.rt.df$Level_Type = 'PHR'
state.rt.df$Level_Type = 'State'
combined.rt.df <- rbind(subset(county.rt.df, Date != max(Date)),
subset(metro.rt.df, Date != max(Date)),
subset(TSA.rt.df, Date != max(Date)),
subset(state.rt.df, Date != max(Date)),
subset(phr.rt.df, Date != max(Date)))
write.csv(combined.rt.df, 'statistical-output/rt/stacked_rt.csv', row.names = FALSE)
write.csv(combined.rt.df, 'tableau/stacked_rt.csv', row.names = FALSE)TIME SERIES
#ts.rt.data.clean function will clean data, dropping unwanted variables
ts.data.clean<-function(covid.data) {
ts_cols <- c("County", "Date", "TSA", "TSA_Name", "PHR", "PHR_Name",
"Metro_Area", "Cases_Daily", "Cases_Cumulative",
"Population_DSHS", "Hospitalizations_Total")
covid.data <- covid.data[, names(covid.data) %in% ts_cols]
# set correct types
covid.data$Date <- as.Date(covid.data$Date)
covid.data$Cases_Daily <- as.numeric(covid.data$Cases_Daily)
#change any Cases_Daily below zero to zero
covid.data <- covid.data %>% mutate(Cases_Daily = ifelse(Cases_Daily < 0, 0, Cases_Daily))
return(covid.data)
}#separate county, metro, TSA and State data into separate Data frames
ts.data.organize<-function(mycounty, mytsa, myphr){
### COUNTY ###
# clean county vals and restrict to first date of collection
cleaned.county <- ts.data.clean(mycounty)
cleaned.county <- subset(cleaned.county, !is.na(Date) & Date >= as.Date('2020-03-04'))
# Set column names
county.df <- cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area)
county.df$Level <- 'County'
colnames(county.df)[2] = 'Level_Name'
#change hospitalizations to numeric
mytsa$Hospitalizations_Total<-as.numeric(mytsa$Hospitalizations_Total)
# PHR
PHR.df <- ts.data.clean(myphr)
PHR.df$Level = 'PHR'
colnames(PHR.df)[3] = 'Level_Name'
### METRO ###
metro.temp<-cleaned.county %>%
group_by(Metro_Area, Date) %>%
mutate(metro_Cases_Daily = sum(Cases_Daily, na.rm = TRUE)) %>%
mutate(metro_Cases_Cumulative = sum(Cases_Cumulative, na.rm=TRUE)) %>%
mutate(metro_pop_DSHS = sum(Population_DSHS, na.rm=TRUE)) %>%
dplyr::select(Date, Metro_Area, metro_Cases_Daily, metro_Cases_Cumulative, metro_pop_DSHS) %>%
distinct()
metro.df <- data.frame(Date = metro.temp$Date,
Level = 'metro',
Level_Name = metro.temp$Metro_Area,
Cases_Daily = metro.temp$metro_Cases_Daily,
Cases_Cumulative = metro.temp$metro_Cases_Cumulative,
Population_DSHS = metro.temp$metro_pop_DSHS)
# drop NA dates
metro.df <- subset(metro.df, !is.na(Date) & Date>=as.Date('2020-03-04'))
### TSA ###
TSA.df <- ts.data.clean(mytsa)
TSA.df <- subset(TSA.df, !is.na(Date) & Date >= as.Date('2020-03-04'))
TSA.df$Level <- 'TSA'
colnames(TSA.df)[2] <- 'Level_Name'
### STATE ###
state.temp<-TSA.df %>%
group_by(Date) %>%
mutate(state_Cases_Daily = sum(Cases_Daily, na.rm = TRUE)) %>%
mutate(state_Cases_Cumulative = sum(Cases_Cumulative, na.rm = TRUE)) %>%
mutate(state_pop_DSHS = sum(Population_DSHS, na.rm = TRUE)) %>%
mutate(state_hosp = sum(Hospitalizations_Total, na.rm=TRUE))%>%
dplyr::select(Date, state_Cases_Daily, state_Cases_Cumulative, state_pop_DSHS, state_hosp) %>%
distinct()
state.df <- data.frame(Date=state.temp$Date,
Level = 'State',
Level_Name = 'Texas',
Cases_Daily=state.temp$state_Cases_Daily,
Cases_Cumulative=state.temp$state_Cases_Cumulative,
Population_DSHS=state.temp$state_pop_DSHS,
Hospitalizations_Total=state.temp$state_hosp)
state.df <- subset(state.df, Date>=as.Date('2020-03-04'))
### OUTPUT ###
ts.df.out <- list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df, phr = PHR.df)
return(ts.df.out)
}# Compute forecast (UPDATE PREDICTION PERIOD [days] AS NEEDED)
covid.arima.forecast<-function(mydata, prediction.period = 10, mindate) {
maxdate <- max(mydata$Date)
# mindate <- as.Date('2020-05-01')
pred_start_label = format(mindate, format = '%m_%d')
mydata = subset(mydata, Date >= mindate)
model.length <- as.numeric(length(mydata$Date) + prediction.period)
case_avg_date = as.Date('2020-05-01')
if(max(mydata$Cases_Daily>=100, na.rm = TRUE)) {
# arima requires cases to be a timeseries vector
# convert daily cases to time series
my.timeseries<-ts(mydata$Cases_Daily)
print(mydata[1,2])
#load package(pracma)
library(pracma)
my.timeseries<-movavg(my.timeseries,7,"s")
#d=0 restricts first differencing to 0 so that daily cases aren't differenced
arima.fit <- forecast::auto.arima(my.timeseries)
# obtain diagnostic plots for ideal arima (p,d,q) selection
acf <- ggAcf(my.timeseries, lag.max=30)
pacf <- ggPacf(my.timeseries, lag.max=30)
ts.diagnostics <- grid.arrange(acf, pacf, nrow=2)
ggsave(paste0('statistical-output/time-series/diagnostics/',
mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'),
plot = ts.diagnostics)
# save parameters from arima autofit
p<- arima.fit$arma[1] # autoregressive order
q<- arima.fit$arma[2] # moving average order
d<-arima.fit$arma[6] # differencing order from model
# 10 day forecast, CI for lower and upper has confidence level 70% set by level =c(70,70)
arima.forecast <- forecast::forecast(arima.fit, h = prediction.period, level=c(70,70))
#return a dataframe of the arima model(Daily cases by date)
arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
Cases_Raw = c(mydata$Cases_Daily, rep(NA, times = prediction.period)),
Cases_Daily = c(my.timeseries, arima.forecast[['mean']]),
CI_Lower = c(rep(NA, times = length(my.timeseries)),
arima.forecast[['lower']][, 2]),
CI_Upper = c(rep(NA, times = length(my.timeseries)),
arima.forecast[['upper']][, 2]))
# Order_AutoReg = c(rep(p, times = model.length)),
# Order_Moving_Avg = c(rep(q, times = model.length)),
# Differencing = c(rep(d, times = model.length)))
# save prediction plot for preliminary review
arima.plot <- ggplot(arima.out, aes(x=Date, y = Cases_Daily))+
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), fill = "red", alpha = 0.5, size = 0.1)+
geom_line(color = "black", size = 1)+
labs(y = 'Daily Cases (7-Day Moving Average)', x = 'Date',
title = mydata$Level_Name[1]) +
scale_x_date(limits = c(mindate, maxdate + prediction.period),
date_labels = '%b-%d', date_breaks = '1 week') +
ggpubr::theme_pubr() +
theme(axis.text.x = element_text(angle = -90))
ggsave(plot = arima.plot, paste0('statistical-output/time-series/plots/',
mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'))
} else {
# insufficient data catch: return NA values for predictions
arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
Cases_Raw = c(mydata$Cases_Daily, rep(NA, times = prediction.period)),
Cases_Daily = rep(NA, times = model.length),
CI_Lower = rep(NA, times = model.length),
CI_Upper = rep(NA, times = model.length))
# Order_AutoReg = c(rep(NA, times = model.length)),
# Order_Moving_Avg = c(rep(NA, times = model.length)),
# Differencing = c(rep(NA, times= model.length)))
}
#replace CI lower limit with 0 if negative
arima.out$CI_Lower <- ifelse(arima.out$CI_Lower>=0 ,arima.out$CI_Lower, 0)
return(arima.out)
}County diangostics
mydata = county.Daily
maxdate = max(mydata$Date)
case_features = subset(mydata, Date >= maxdate-14) %>%
group_by(Level_Name) %>%
mutate(criteria1 = mean(Cases_Daily)) %>%
mutate(criteria2 = mean(Cases_Cumulative)) %>%
mutate(criteria3 = (mean(Cases_Cumulative) / Population_DSHS)) %>%
dplyr::select(Level_Name, Date, criteria1, criteria2, criteria3, PHR, TSA)
criteria1_vals = subset(case_features, Date == maxdate) %>%
group_by(Date) %>%
do(data.frame(t(quantile(.$criteria1, probs = c(0.4, 0.5, 0.6, 0.7, 0.8)))))
criteria1_features = merge(criteria1_vals, case_features, by = 'Date')
ggplot(subset(criteria1_features, Date == maxdate), aes(y = Level_Name, x = criteria1)) +
geom_point() +
geom_vline(xintercept = criteria1_features$X40[1], color = 'red') +
geom_vline(xintercept = criteria1_features$X50.[1], color = 'orange') +
geom_vline(xintercept = criteria1_features$X60.[1], color = 'green') +
geom_vline(xintercept = criteria1_features$X70.[1], color = 'blue') +
geom_vline(xintercept = criteria1_features$X80.[1], color = 'purple') +
labs(title = paste0('mean(Cases_Daily) by county quantile (Data from: ', maxdate-14, ')'),
subtitle = paste0('40%: ', round(criteria1_features$X40.[1], 2), ' | ',
'50%: ', round(criteria1_features$X50.[1], 2), ' | ',
'60%: ', round(criteria1_features$X60.[1], 2), ' |',
'70%: ', round(criteria1_features$X70.[1], 2), ' | ',
'80%: ', round(criteria1_features$X80.[1], 2)),
x = '', y ='') +
facet_wrap(~TSA, scales = 'free') +
ggpubr::theme_pubr(border = TRUE) +
theme(axis.text.y = element_text(size = 5)) +
theme(axis.text.x = element_text(size = 6))
ggsave('statistical-output/time-series/diagnostics/threshold1.png', width = 10, height = 10, dpi = 600)
# TODO: add other criteriaCounty
# apply arima across all counties
county.arima.output.1 <- nlme::gapply(county.Daily,
FUN = covid.arima.forecast,
groups = county.Daily$Level_Name,
mindate = as.Date('2020-03-04'))## [1] Anderson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala##
## Attaching package: 'pracma'## The following object is masked from 'package:mgcv':
##
## magic## The following object is masked from 'package:purrr':
##
## cross## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Angelina
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Bell
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Bexar
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Brazoria
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Brazos
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Cameron
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Chambers
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Collin
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Comal
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Dallas
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Denton
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] DeWitt
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Ector
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] El Paso
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Ellis
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Fort Bend
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Galveston
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Gregg
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Grimes
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Harris
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Hays
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Hidalgo
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Jefferson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Johnson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Jones
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Kaufman
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] La Salle
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Lubbock
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Maverick
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] McLennan
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Medina
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Midland
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Montgomery
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Moore
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Nueces
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Orange
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Parker
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Polk
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Potter
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Randall
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Rusk
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Scurry
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Smith
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Starr
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Tarrant
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Titus
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Tom Green
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Travis
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Val Verde
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Victoria
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Walker
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Webb
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Williamson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Wilson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala## Saving 7 x 5 in image
## Saving 7 x 5 in image
# bind list of dataframes to one dataframe
county.arima.df.1 <- data.table::rbindlist(county.arima.output.1, idcol = TRUE)
colnames(county.arima.df.1)[1] <- 'County'
# county.arima.output.2 <- nlme::gapply(county.Daily,
# FUN = covid.arima.forecast,
# groups = county.Daily$Level_Name,
# mindate = as.Date('2020-06-03'))
#
#
# # bind list of dataframes to one dataframe
# county.arima.df.2 <- data.table::rbindlist(county.arima.output.2, idcol = TRUE)
# colnames(county.arima.df.2)[1] <- 'County'Metro
# apply arima across both metro values
metro.arima.output.1 <- nlme::gapply(metro.Daily,
FUN = covid.arima.forecast,
groups = metro.Daily$Level_Name,
mindate = as.Date('2020-03-04'))## [1] metro
## Levels: metro## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] metro
## Levels: metro## Saving 7 x 5 in image
## Saving 7 x 5 in image
# bind list of dataframes to one dataframe
metro.arima.df.1 <- data.table::rbindlist(metro.arima.output.1, idcol = TRUE)
colnames(metro.arima.df.1)[1] <- 'Metro_Area'
#
# metro.arima.output.2 <- nlme::gapply(metro.Daily,
# FUN = covid.arima.forecast,
# groups = metro.Daily$Level_Name,
# mindate = as.Date('2020-06-03'))
#
# # bind list of dataframes to one dataframe
# metro.arima.df.2 <- data.table::rbindlist(metro.arima.output.2, idcol = TRUE)
# colnames(metro.arima.df.2)[1] <- 'Metro_Area'TSA
# All data
TSA.arima.output.1 <- nlme::gapply(TSA.Daily,
FUN = covid.arima.forecast,
groups = TSA.Daily$Level_Name,
mindate = as.Date('2020-03-04'))## [1] A
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] B
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] D
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] E
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] F
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] G
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] H
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] I
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] J
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] K
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] L
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] M
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] N
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] O
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] P
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] Q
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] R
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] S
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] T
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] U
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] V
## Levels: A B C D E F G H I J K L M N O P Q R S T U V## Saving 7 x 5 in image
## Saving 7 x 5 in image
TSA.arima.df.1 <- data.table::rbindlist(TSA.arima.output.1, idcol=TRUE)
colnames(TSA.arima.df.1)[1] <- 'TSA'
TSA.arima.df.1 <- merge(TSA.arima.df.1, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))
# combine TSA
TSA.arima.df.1$TSA <- paste0(TSA.arima.df.1$TSA, ' - ', TSA.arima.df.1$TSA_Name)
TSA.arima.df.1$TSA_Name <- NULLPHR
PHR.arima.output <- nlme::gapply(PHR.Daily,
FUN = covid.arima.forecast,
groups = PHR.Daily$Level_Name,
mindate = as.Date('2020-03-04'))## [1] 2/3
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 9/10
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 11
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 6/5S
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 1
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 8
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 7
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
## [1] 4/5N
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10## Saving 7 x 5 in image
## Saving 7 x 5 in image
PHR.arima.df <- data.table::rbindlist(PHR.arima.output, idcol=TRUE)
colnames(PHR.arima.df)[1] <- 'Level_Name'
PHR.arima.df <- merge(PHR.arima.df, unique(PHR.Daily[, c('PHR', 'Level_Name')]), by = c('Level_Name'))
# combine PHR
PHR.arima.df$Level_Name <- paste0(PHR.arima.df$PHR, ' - ', PHR.arima.df$Level_Name)
PHR.arima.df$PHR <- NULL
colnames(PHR.arima.df)[1] <- 'PHR'State
texas.arima.df.1 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-03-04'))## Saving 7 x 5 in image
## Saving 7 x 5 in image
# texas.arima.df.2 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-06-03'))#Create plot function, forecast.data is the dataframe containing true data
# and the forecast data. days.forecast is the number of days that are forecast
# in the dataframe (here we are forecasting 10 days, so it should be 10 unless we
# change this at some point down the road)
forecast.plot<-function(forecast.data, plot.title, y.label){
mindate<-as.POSIXct("2020-03-04")
maxdate<-as.POSIXct(max(forecast.data$Date))
ts.plot<-ggplot(aes(x=as.POSIXct(Date), y=Cases_Daily), data=forecast.data)+
geom_ribbon(aes(ymin=CI_Lower, ymax=CI_Upper), fill="grey50", size=0.1)+
geom_line(color="blue", size=1)+
geom_line(aes(x=as.POSIXct(Date), y=CI_Lower), color="grey", size=0.1)+
geom_line(aes(x=as.POSIXct(Date),y=CI_Upper), color="grey", size=0.1)+
scale_x_datetime(limits = c(mindate, maxdate))+
xlab("Date")+
ylab(y.label)+
#Can use the following title if we are running using nlme for data frame
#ggtitle(paste(forecast.data$.id,": TS Daily Cases + 10 Day Forcast",sep=""))
ggtitle(plot.title)
ts.plot
}######### View Output Table & Graph for TSA Q ##########
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows#subset TSA Q arima data
TSA_Q.arima.df<-TSA.arima.df.1%>%subset(TSA.arima.df.1$TSA== "Q - Houston")
#apply the plot function to TSA Q arima data frame
TSA_Q.arima.df%>%kable(caption="Greater Houston ARIMA - Daily Cases")%>%kable_styling(full_width=FALSE)| TSA | Date | Cases_Raw | Cases_Daily | CI_Lower | CI_Upper |
|---|---|---|---|---|---|
| Q - Houston | 2020-03-04 | 0 | 0.000000 | NA | NA |
| Q - Houston | 2020-03-05 | 0 | 0.000000 | NA | NA |
| Q - Houston | 2020-03-06 | 5 | 1.666667 | NA | NA |
| Q - Houston | 2020-03-07 | NA | NA | NA | NA |
| Q - Houston | 2020-03-08 | NA | NA | NA | NA |
| Q - Houston | 2020-03-09 | 6 | NA | NA | NA |
| Q - Houston | 2020-03-10 | 0 | NA | NA | NA |
| Q - Houston | 2020-03-11 | 0 | NA | NA | NA |
| Q - Houston | 2020-03-12 | 3 | NA | NA | NA |
| Q - Houston | 2020-03-13 | 0 | NA | NA | NA |
| Q - Houston | 2020-03-14 | NA | NA | NA | NA |
| Q - Houston | 2020-03-15 | 9 | NA | NA | NA |
| Q - Houston | 2020-03-16 | 0 | NA | NA | NA |
| Q - Houston | 2020-03-17 | 0 | NA | NA | NA |
| Q - Houston | 2020-03-18 | 0 | NA | NA | NA |
| Q - Houston | 2020-03-19 | 2 | NA | NA | NA |
| Q - Houston | 2020-03-20 | 17 | NA | NA | NA |
| Q - Houston | 2020-03-21 | 2 | 4.285714 | NA | NA |
| Q - Houston | 2020-03-22 | 1 | 3.142857 | NA | NA |
| Q - Houston | 2020-03-23 | 0 | 3.142857 | NA | NA |
| Q - Houston | 2020-03-24 | 60 | 11.714286 | NA | NA |
| Q - Houston | 2020-03-25 | 109 | 27.285714 | NA | NA |
| Q - Houston | 2020-03-26 | 72 | 37.285714 | NA | NA |
| Q - Houston | 2020-03-27 | 60 | 43.428571 | NA | NA |
| Q - Houston | 2020-03-28 | 48 | 50.000000 | NA | NA |
| Q - Houston | 2020-03-29 | 242 | 84.428571 | NA | NA |
| Q - Houston | 2020-03-30 | 102 | 99.000000 | NA | NA |
| Q - Houston | 2020-03-31 | 52 | 97.857143 | NA | NA |
| Q - Houston | 2020-04-01 | 190 | 109.428571 | NA | NA |
| Q - Houston | 2020-04-02 | 227 | 131.571429 | NA | NA |
| Q - Houston | 2020-04-03 | 138 | 142.714286 | NA | NA |
| Q - Houston | 2020-04-04 | 212 | 166.142857 | NA | NA |
| Q - Houston | 2020-04-05 | 204 | 160.714286 | NA | NA |
| Q - Houston | 2020-04-06 | 138 | 165.857143 | NA | NA |
| Q - Houston | 2020-04-07 | 473 | 226.000000 | NA | NA |
| Q - Houston | 2020-04-08 | 477 | 267.000000 | NA | NA |
| Q - Houston | 2020-04-09 | 231 | 267.571429 | NA | NA |
| Q - Houston | 2020-04-10 | 781 | 359.428571 | NA | NA |
| Q - Houston | 2020-04-11 | 275 | 368.428571 | NA | NA |
| Q - Houston | 2020-04-12 | 333 | 386.857143 | NA | NA |
| Q - Houston | 2020-04-13 | 76 | 378.000000 | NA | NA |
| Q - Houston | 2020-04-14 | 163 | 333.714286 | NA | NA |
| Q - Houston | 2020-04-15 | 292 | 307.285714 | NA | NA |
| Q - Houston | 2020-04-16 | 236 | 308.000000 | NA | NA |
| Q - Houston | 2020-04-17 | 267 | 234.571429 | NA | NA |
| Q - Houston | 2020-04-18 | 277 | 234.857143 | NA | NA |
| Q - Houston | 2020-04-19 | 211 | 217.428571 | NA | NA |
| Q - Houston | 2020-04-20 | 194 | 234.285714 | NA | NA |
| Q - Houston | 2020-04-21 | 216 | 241.857143 | NA | NA |
| Q - Houston | 2020-04-22 | 201 | 228.857143 | NA | NA |
| Q - Houston | 2020-04-23 | 195 | 223.000000 | NA | NA |
| Q - Houston | 2020-04-24 | 181 | 210.714286 | NA | NA |
| Q - Houston | 2020-04-25 | 214 | 201.714286 | NA | NA |
| Q - Houston | 2020-04-26 | 171 | 196.000000 | NA | NA |
| Q - Houston | 2020-04-27 | 161 | 191.285714 | NA | NA |
| Q - Houston | 2020-04-28 | 139 | 180.285714 | NA | NA |
| Q - Houston | 2020-04-29 | 219 | 182.857143 | NA | NA |
| Q - Houston | 2020-04-30 | 276 | 194.428571 | NA | NA |
| Q - Houston | 2020-05-01 | 275 | 207.857143 | NA | NA |
| Q - Houston | 2020-05-02 | 251 | 213.142857 | NA | NA |
| Q - Houston | 2020-05-03 | 243 | 223.428571 | NA | NA |
| Q - Houston | 2020-05-04 | 182 | 226.428571 | NA | NA |
| Q - Houston | 2020-05-05 | 136 | 226.000000 | NA | NA |
| Q - Houston | 2020-05-06 | 257 | 231.428571 | NA | NA |
| Q - Houston | 2020-05-07 | 181 | 217.857143 | NA | NA |
| Q - Houston | 2020-05-08 | 213 | 209.000000 | NA | NA |
| Q - Houston | 2020-05-09 | 263 | 210.714286 | NA | NA |
| Q - Houston | 2020-05-10 | 224 | 208.000000 | NA | NA |
| Q - Houston | 2020-05-11 | 92 | 195.142857 | NA | NA |
| Q - Houston | 2020-05-12 | 344 | 224.857143 | NA | NA |
| Q - Houston | 2020-05-13 | 310 | 232.428571 | NA | NA |
| Q - Houston | 2020-05-14 | 240 | 240.857143 | NA | NA |
| Q - Houston | 2020-05-15 | 243 | 245.142857 | NA | NA |
| Q - Houston | 2020-05-16 | 332 | 255.000000 | NA | NA |
| Q - Houston | 2020-05-17 | 141 | 243.142857 | NA | NA |
| Q - Houston | 2020-05-18 | 343 | 279.000000 | NA | NA |
| Q - Houston | 2020-05-19 | 214 | 260.428571 | NA | NA |
| Q - Houston | 2020-05-20 | 391 | 272.000000 | NA | NA |
| Q - Houston | 2020-05-21 | 167 | 261.571429 | NA | NA |
| Q - Houston | 2020-05-22 | 243 | 261.571429 | NA | NA |
| Q - Houston | 2020-05-23 | 288 | 255.285714 | NA | NA |
| Q - Houston | 2020-05-24 | 275 | 274.428571 | NA | NA |
| Q - Houston | 2020-05-25 | 171 | 249.857143 | NA | NA |
| Q - Houston | 2020-05-26 | 126 | 237.285714 | NA | NA |
| Q - Houston | 2020-05-27 | 561 | 261.571429 | NA | NA |
| Q - Houston | 2020-05-28 | 463 | 303.857143 | NA | NA |
| Q - Houston | 2020-05-29 | 262 | 306.571429 | NA | NA |
| Q - Houston | 2020-05-30 | 379 | 319.571429 | NA | NA |
| Q - Houston | 2020-05-31 | 752 | 387.714286 | NA | NA |
| Q - Houston | 2020-06-01 | 83 | 375.142857 | NA | NA |
| Q - Houston | 2020-06-02 | 526 | 432.285714 | NA | NA |
| Q - Houston | 2020-06-03 | 603 | 438.285714 | NA | NA |
| Q - Houston | 2020-06-04 | 375 | 425.714286 | NA | NA |
| Q - Houston | 2020-06-05 | 484 | 457.428571 | NA | NA |
| Q - Houston | 2020-06-06 | 446 | 467.000000 | NA | NA |
| Q - Houston | 2020-06-07 | 489 | 429.428571 | NA | NA |
| Q - Houston | 2020-06-08 | 163 | 440.857143 | NA | NA |
| Q - Houston | 2020-06-09 | 416 | 425.142857 | NA | NA |
| Q - Houston | 2020-06-10 | 424 | 399.571429 | NA | NA |
| Q - Houston | 2020-06-11 | 455 | 411.000000 | NA | NA |
| Q - Houston | 2020-06-12 | 385 | 396.857143 | NA | NA |
| Q - Houston | 2020-06-13 | 427 | 394.142857 | NA | NA |
| Q - Houston | 2020-06-14 | 415 | 383.571429 | NA | NA |
| Q - Houston | 2020-06-15 | 221 | 391.857143 | NA | NA |
| Q - Houston | 2020-06-16 | 572 | 414.142857 | NA | NA |
| Q - Houston | 2020-06-17 | 606 | 440.142857 | NA | NA |
| Q - Houston | 2020-06-18 | 708 | 476.285714 | NA | NA |
| Q - Houston | 2020-06-19 | 544 | 499.000000 | NA | NA |
| Q - Houston | 2020-06-20 | 1443 | 644.142857 | NA | NA |
| Q - Houston | 2020-06-21 | 1224 | 759.714286 | NA | NA |
| Q - Houston | 2020-06-22 | 305 | 771.714286 | NA | NA |
| Q - Houston | 2020-06-23 | 2187 | 1002.428571 | NA | NA |
| Q - Houston | 2020-06-24 | 1581 | 1141.714286 | NA | NA |
| Q - Houston | 2020-06-25 | 1575 | 1265.571429 | NA | NA |
| Q - Houston | 2020-06-26 | 1438 | 1393.285714 | NA | NA |
| Q - Houston | 2020-06-27 | 1604 | 1416.285714 | NA | NA |
| Q - Houston | 2020-06-28 | 1019 | 1387.000000 | NA | NA |
| Q - Houston | 2020-06-29 | 142 | 1363.714286 | NA | NA |
| Q - Houston | 2020-06-30 | 1569 | 1275.428571 | NA | NA |
| Q - Houston | 2020-07-01 | 953 | 1185.714286 | NA | NA |
| Q - Houston | 2020-07-02 | 1713 | 1205.428571 | NA | NA |
| Q - Houston | 2020-07-03 | 1554 | 1222.000000 | NA | NA |
| Q - Houston | 2020-07-04 | 1356 | 1186.571429 | NA | NA |
| Q - Houston | 2020-07-05 | 594 | 1125.857143 | NA | NA |
| Q - Houston | 2020-07-06 | 734 | 1210.428571 | NA | NA |
| Q - Houston | 2020-07-07 | 1580 | 1212.000000 | NA | NA |
| Q - Houston | 2020-07-08 | 1785 | 1330.857143 | NA | NA |
| Q - Houston | 2020-07-09 | 996 | 1228.428571 | NA | NA |
| Q - Houston | 2020-07-10 | 1250 | 1185.000000 | NA | NA |
| Q - Houston | 2020-07-11 | 1383 | 1188.857143 | NA | NA |
| Q - Houston | 2020-07-12 | 2137 | 1409.285714 | NA | NA |
| Q - Houston | 2020-07-13 | 1524 | 1522.142857 | NA | NA |
| Q - Houston | 2020-07-14 | 2595 | 1667.142857 | NA | NA |
| Q - Houston | 2020-07-15 | 2253 | 1734.000000 | NA | NA |
| Q - Houston | 2020-07-16 | 2314 | 1922.285714 | NA | NA |
| Q - Houston | 2020-07-17 | 1935 | 2020.142857 | NA | NA |
| Q - Houston | 2020-07-18 | 2205 | 2137.571429 | NA | NA |
| Q - Houston | 2020-07-19 | 1584 | 2058.571429 | NA | NA |
| Q - Houston | 2020-07-20 | 1031 | 1988.142857 | NA | NA |
| Q - Houston | 2020-07-21 | 1780 | 1871.714286 | NA | NA |
| Q - Houston | 2020-07-22 | 1685 | 1790.571429 | NA | NA |
| Q - Houston | 2020-07-23 | 1725 | 1706.428571 | NA | NA |
| Q - Houston | 2020-07-24 | NA | 1685.220532 | 1635.552 | 1734.889 |
| Q - Houston | 2020-07-25 | NA | 1644.817533 | 1556.960 | 1732.675 |
| Q - Houston | 2020-07-26 | NA | 1645.979416 | 1518.047 | 1773.911 |
| Q - Houston | 2020-07-27 | NA | 1622.907440 | 1461.045 | 1784.770 |
| Q - Houston | 2020-07-28 | NA | 1631.026703 | 1436.057 | 1825.997 |
| Q - Houston | 2020-07-29 | NA | 1615.438692 | 1392.033 | 1838.844 |
| Q - Houston | 2020-07-30 | NA | 1625.044219 | 1374.095 | 1875.993 |
| Q - Houston | 2020-07-31 | NA | 1613.179261 | 1338.034 | 1888.324 |
| Q - Houston | 2020-08-01 | NA | 1622.414821 | 1323.723 | 1921.107 |
| Q - Houston | 2020-08-02 | NA | 1612.760978 | 1292.945 | 1932.577 |
#nlme::gapply(TSA.arima.df, FUN = forecast.plot, groups=TSA.daily$Level_Name)
forecast.plot(TSA_Q.arima.df, "TSA Q - Greater Houston Daily Cases", "Daily Cases")## Warning: Removed 1 rows containing missing values (geom_path).## Warning: Removed 142 rows containing missing values (geom_path).
## Warning: Removed 142 rows containing missing values (geom_path).
######### View Output Table & Graph for Texas ##########
texas.arima.df.1%>%kable(caption="Texas Arima - Daily Cases")%>%kable_styling(full_width=FALSE)| Date | Cases_Raw | Cases_Daily | CI_Lower | CI_Upper |
|---|---|---|---|---|
| 2020-03-04 | 0 | 0.000000 | NA | NA |
| 2020-03-05 | 0 | 0.000000 | NA | NA |
| 2020-03-06 | 5 | 1.666667 | NA | NA |
| 2020-03-07 | 0 | 1.250000 | NA | NA |
| 2020-03-08 | 0 | 1.000000 | NA | NA |
| 2020-03-09 | 7 | 2.000000 | NA | NA |
| 2020-03-10 | 3 | 2.142857 | NA | NA |
| 2020-03-11 | 3 | 2.571429 | NA | NA |
| 2020-03-12 | 4 | 3.142857 | NA | NA |
| 2020-03-13 | 0 | 2.428571 | NA | NA |
| 2020-03-14 | 0 | 2.428571 | NA | NA |
| 2020-03-15 | 34 | 7.285714 | NA | NA |
| 2020-03-16 | 0 | 6.285714 | NA | NA |
| 2020-03-17 | 7 | 6.857143 | NA | NA |
| 2020-03-18 | 19 | 9.142857 | NA | NA |
| 2020-03-19 | 26 | 12.285714 | NA | NA |
| 2020-03-20 | 67 | 21.857143 | NA | NA |
| 2020-03-21 | 60 | 30.428571 | NA | NA |
| 2020-03-22 | 28 | 29.571429 | NA | NA |
| 2020-03-23 | 24 | 33.000000 | NA | NA |
| 2020-03-24 | 425 | 92.714286 | NA | NA |
| 2020-03-25 | 263 | 127.571429 | NA | NA |
| 2020-03-26 | 419 | 183.714286 | NA | NA |
| 2020-03-27 | 337 | 222.285714 | NA | NA |
| 2020-03-28 | 317 | 259.000000 | NA | NA |
| 2020-03-29 | 504 | 327.000000 | NA | NA |
| 2020-03-30 | 322 | 369.571429 | NA | NA |
| 2020-03-31 | 392 | 364.857143 | NA | NA |
| 2020-04-01 | 730 | 431.571429 | NA | NA |
| 2020-04-02 | 669 | 467.285714 | NA | NA |
| 2020-04-03 | 659 | 513.285714 | NA | NA |
| 2020-04-04 | 788 | 580.571429 | NA | NA |
| 2020-04-05 | 681 | 605.857143 | NA | NA |
| 2020-04-06 | 484 | 629.000000 | NA | NA |
| 2020-04-07 | 988 | 714.142857 | NA | NA |
| 2020-04-08 | 1092 | 765.857143 | NA | NA |
| 2020-04-09 | 877 | 795.571429 | NA | NA |
| 2020-04-10 | 1441 | 907.285714 | NA | NA |
| 2020-04-11 | 890 | 921.857143 | NA | NA |
| 2020-04-12 | 923 | 956.428571 | NA | NA |
| 2020-04-13 | 422 | 947.571429 | NA | NA |
| 2020-04-14 | 718 | 909.000000 | NA | NA |
| 2020-04-15 | 868 | 877.000000 | NA | NA |
| 2020-04-16 | 963 | 889.285714 | NA | NA |
| 2020-04-17 | 916 | 814.285714 | NA | NA |
| 2020-04-18 | 889 | 814.142857 | NA | NA |
| 2020-04-19 | 663 | 777.000000 | NA | NA |
| 2020-04-20 | 535 | 793.142857 | NA | NA |
| 2020-04-21 | 738 | 796.000000 | NA | NA |
| 2020-04-22 | 873 | 796.714286 | NA | NA |
| 2020-04-23 | 876 | 784.285714 | NA | NA |
| 2020-04-24 | 862 | 776.571429 | NA | NA |
| 2020-04-25 | 968 | 787.857143 | NA | NA |
| 2020-04-26 | 859 | 815.857143 | NA | NA |
| 2020-04-27 | 666 | 834.571429 | NA | NA |
| 2020-04-28 | 874 | 854.000000 | NA | NA |
| 2020-04-29 | 885 | 855.714286 | NA | NA |
| 2020-04-30 | 1033 | 878.142857 | NA | NA |
| 2020-05-01 | 1142 | 918.142857 | NA | NA |
| 2020-05-02 | 1293 | 964.571429 | NA | NA |
| 2020-05-03 | 1026 | 988.428571 | NA | NA |
| 2020-05-04 | 784 | 1005.285714 | NA | NA |
| 2020-05-05 | 1037 | 1028.571429 | NA | NA |
| 2020-05-06 | 1053 | 1052.571429 | NA | NA |
| 2020-05-07 | 1135 | 1067.142857 | NA | NA |
| 2020-05-08 | 1219 | 1078.142857 | NA | NA |
| 2020-05-09 | 1251 | 1072.142857 | NA | NA |
| 2020-05-10 | 1009 | 1069.714286 | NA | NA |
| 2020-05-11 | 1000 | 1100.571429 | NA | NA |
| 2020-05-12 | 1179 | 1120.857143 | NA | NA |
| 2020-05-13 | 1355 | 1164.000000 | NA | NA |
| 2020-05-14 | 1448 | 1208.714286 | NA | NA |
| 2020-05-15 | 1347 | 1227.000000 | NA | NA |
| 2020-05-16 | 1801 | 1305.571429 | NA | NA |
| 2020-05-17 | 785 | 1273.571429 | NA | NA |
| 2020-05-18 | 909 | 1260.571429 | NA | NA |
| 2020-05-19 | 1219 | 1266.285714 | NA | NA |
| 2020-05-20 | 1411 | 1274.285714 | NA | NA |
| 2020-05-21 | 945 | 1202.428571 | NA | NA |
| 2020-05-22 | 1241 | 1187.285714 | NA | NA |
| 2020-05-23 | 1060 | 1081.428571 | NA | NA |
| 2020-05-24 | 839 | 1089.142857 | NA | NA |
| 2020-05-25 | 623 | 1048.285714 | NA | NA |
| 2020-05-26 | 620 | 962.714286 | NA | NA |
| 2020-05-27 | 1361 | 955.571429 | NA | NA |
| 2020-05-28 | 1855 | 1085.571429 | NA | NA |
| 2020-05-29 | 1230 | 1084.000000 | NA | NA |
| 2020-05-30 | 1333 | 1123.000000 | NA | NA |
| 2020-05-31 | 1949 | 1281.571429 | NA | NA |
| 2020-06-01 | 593 | 1277.285714 | NA | NA |
| 2020-06-02 | 1688 | 1429.857143 | NA | NA |
| 2020-06-03 | 1703 | 1478.714286 | NA | NA |
| 2020-06-04 | 1650 | 1449.428571 | NA | NA |
| 2020-06-05 | 1693 | 1515.571429 | NA | NA |
| 2020-06-06 | 1942 | 1602.571429 | NA | NA |
| 2020-06-07 | 1426 | 1527.857143 | NA | NA |
| 2020-06-08 | 638 | 1534.285714 | NA | NA |
| 2020-06-09 | 1853 | 1557.857143 | NA | NA |
| 2020-06-10 | 2509 | 1673.000000 | NA | NA |
| 2020-06-11 | 1826 | 1698.142857 | NA | NA |
| 2020-06-12 | 2097 | 1755.857143 | NA | NA |
| 2020-06-13 | 2331 | 1811.428571 | NA | NA |
| 2020-06-14 | 1843 | 1871.000000 | NA | NA |
| 2020-06-15 | 1254 | 1959.000000 | NA | NA |
| 2020-06-16 | 4098 | 2279.714286 | NA | NA |
| 2020-06-17 | 3140 | 2369.857143 | NA | NA |
| 2020-06-18 | 3516 | 2611.285714 | NA | NA |
| 2020-06-19 | 3454 | 2805.142857 | NA | NA |
| 2020-06-20 | 4430 | 3105.000000 | NA | NA |
| 2020-06-21 | 3866 | 3394.000000 | NA | NA |
| 2020-06-22 | 3280 | 3683.428571 | NA | NA |
| 2020-06-23 | 5522 | 3886.857143 | NA | NA |
| 2020-06-24 | 5551 | 4231.285714 | NA | NA |
| 2020-06-25 | 5996 | 4585.571429 | NA | NA |
| 2020-06-26 | 5707 | 4907.428571 | NA | NA |
| 2020-06-27 | 5742 | 5094.857143 | NA | NA |
| 2020-06-28 | 5357 | 5307.857143 | NA | NA |
| 2020-06-29 | 4288 | 5451.857143 | NA | NA |
| 2020-06-30 | 6975 | 5659.428571 | NA | NA |
| 2020-07-01 | 8076 | 6020.142857 | NA | NA |
| 2020-07-02 | 7915 | 6294.285714 | NA | NA |
| 2020-07-03 | 7555 | 6558.285714 | NA | NA |
| 2020-07-04 | 8258 | 6917.714286 | NA | NA |
| 2020-07-05 | 3449 | 6645.142857 | NA | NA |
| 2020-07-06 | 5319 | 6792.428571 | NA | NA |
| 2020-07-07 | 10028 | 7228.571429 | NA | NA |
| 2020-07-08 | 9979 | 7500.428571 | NA | NA |
| 2020-07-09 | 9782 | 7767.142857 | NA | NA |
| 2020-07-10 | 9765 | 8082.857143 | NA | NA |
| 2020-07-11 | 10351 | 8381.857143 | NA | NA |
| 2020-07-12 | 8196 | 9060.000000 | NA | NA |
| 2020-07-13 | 5655 | 9108.000000 | NA | NA |
| 2020-07-14 | 10745 | 9210.428571 | NA | NA |
| 2020-07-15 | 9550 | 9149.142857 | NA | NA |
| 2020-07-16 | 10291 | 9221.857143 | NA | NA |
| 2020-07-17 | 14916 | 9957.714286 | NA | NA |
| 2020-07-18 | 10344 | 9956.714286 | NA | NA |
| 2020-07-19 | 7322 | 9831.857143 | NA | NA |
| 2020-07-20 | 7404 | 10081.714286 | NA | NA |
| 2020-07-21 | 9389 | 9888.000000 | NA | NA |
| 2020-07-22 | 9879 | 9935.000000 | NA | NA |
| 2020-07-23 | 9507 | 9823.000000 | NA | NA |
| 2020-07-24 | NA | 9872.470181 | 9750.871 | 9994.069 |
| 2020-07-25 | NA | 9921.940362 | 9732.133 | 10111.748 |
| 2020-07-26 | NA | 9971.410544 | 9716.680 | 10226.141 |
| 2020-07-27 | NA | 10020.880725 | 9700.671 | 10341.090 |
| 2020-07-28 | NA | 10070.350906 | 9682.882 | 10457.820 |
| 2020-07-29 | NA | 10119.821087 | 9662.823 | 10576.819 |
| 2020-07-30 | NA | 10169.291268 | 9640.287 | 10698.296 |
| 2020-07-31 | NA | 10218.761449 | 9615.194 | 10822.329 |
| 2020-08-01 | NA | 10268.231631 | 9587.530 | 10948.934 |
| 2020-08-02 | NA | 10317.701812 | 9557.309 | 11078.095 |
forecast.plot(texas.arima.df.1, "Texas Daily Cases", "Daily Cases")## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 142 rows containing missing values (geom_path).
## Warning: Removed 142 rows containing missing values (geom_path).
export
write.csv(texas.arima.df.1, 'statistical-output/time-series/state_case_timeseries.csv', row.names = FALSE)
write.csv(TSA.arima.df.1, 'statistical-output/time-series/tsa_case_timeseries.csv', row.names = FALSE)
write.csv(county.arima.df.1, 'statistical-output/time-series/county_case_timeseries.csv', row.names = FALSE)
write.csv(metro.arima.df.1, 'statistical-output/time-series/metro_case_timeseries.csv', row.names = FALSE)
write.csv(PHR.arima.df, 'statistical-output/time-series/phr_case_timeseries.csv', row.names = FALSE)
# write.csv(texas.arima.df.2, 'statistical-output/time-series/state_timeseries_06_03.csv', row.names = FALSE)
# write.csv(TSA.arima.df.2, 'statistical-output/time-series/tsa_timeseries_06_03.csv', row.names = FALSE)
# write.csv(county.arima.df.2, 'statistical-output/time-series/county_timeseries_06_03.csv', row.names = FALSE)
# write.csv(metro.arima.df.2, 'statistical-output/time-series/metro_timeseries_06_03.csv', row.names = FALSE)grouping
colnames(county.arima.df.1)[1] <- 'Level'
colnames(metro.arima.df.1)[1] <- 'Level'
colnames(TSA.arima.df.1)[1] <- 'Level'
texas.arima.df.1$Level <- 'Texas'
colnames(PHR.arima.df)[1] <- 'Level'
county.arima.df.1$Level_Type = 'County'
metro.arima.df.1$Level_Type = 'Metro'
TSA.arima.df.1$Level_Type = 'TSA'
PHR.arima.df$Level_Type = 'PHR'
texas.arima.df.1$Level_Type = 'State'
combined.arima.df.1 <- rbind(county.arima.df.1, metro.arima.df.1, TSA.arima.df.1, texas.arima.df.1, PHR.arima.df)
write.csv(combined.arima.df.1, 'statistical-output/time-series/stacked_case_timeseries.csv',
row.names = FALSE)
write.csv(combined.arima.df.1, 'tableau/stacked_case_timeseries.csv',
row.names = FALSE)
#
# combined.arima.df.2 <- rbind(county.arima.df.2, metro.arima.df.2, TSA.arima.df.2, texas.arima.df.2)
# write.csv(combined.arima.df.2, 'statistical-output/time-series/stacked_timeseries_06_03.csv',
# row.names = FALSE)
#
# write.csv(combined.arima.df.2, 'tableau/stacked_timeseries_06_03.csv',
# row.names = FALSE)Hospitalization Time Series
# Compute forecast (UPDATE PREDICTION PERIOD [days] AS NEEDED)
covid.arima.forecast<-function(mydata, prediction.period = 10, mindate)
{
maxdate <- max(mydata$Date)
# mindate <- as.Date('2020-05-01')
pred_start_label = format(mindate, format = '%m_%d')
mydata = subset(mydata, Date >= mindate)
model.length <- as.numeric(length(mydata$Date) + prediction.period)
if(max(mydata$Hospitalizations_Total>=100, na.rm = TRUE))
{
# arima requires cases to be a timeseries vector
#convert daily cases to time series
my.timeseries<-ts(mydata$Hospitalizations_Total)
#my.timeseries <- rollmeanr(my.timeseries, k=7, na.pad=TRUE, align = 'right')
#load package(pracma)
library(pracma)
my.timeseries<-movavg(my.timeseries,7,"s")
#d=0 restricts first differencing to 0 so that daily cases aren't differenced
arima.fit <- forecast::auto.arima(my.timeseries)
# obtain diagnostic plots for ideal arima (p,d,q) selection
acf <- ggAcf(my.timeseries, lag.max=30)
pacf <- ggPacf(my.timeseries, lag.max=30)
ts.diagnostics <- grid.arrange(acf, pacf, nrow=2)
ggsave(paste0('statistical-output/time-series/diagnostics/',
mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'),
plot = ts.diagnostics)
# save parameters from arima autofit
p<- arima.fit$arma[1] # autoregressive order
q<- arima.fit$arma[2] # moving average order
d<-arima.fit$arma[6] # differencing order from model
# 10 day forecast, CI for lower and upper has confidence level 70% set by level =c(70,70)
arima.forecast <- forecast::forecast(arima.fit, h = prediction.period, level=c(70,70))
#return a dataframe of the arima model(Daily cases by date)
arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
# Cases_Raw = c(mydata$Hospitalizations_Total, rep(NA, times = prediction.period)),
Hospitalizations_Total = c(my.timeseries, arima.forecast[['mean']]),
CI_Lower = c(rep(NA, times = length(my.timeseries)),
arima.forecast[['lower']][, 2]),
CI_Upper = c(rep(NA, times = length(my.timeseries)),
arima.forecast[['upper']][, 2]))
# Order_AutoReg = c(rep(p, times = model.length)),
# Order_Moving_Avg = c(rep(q, times = model.length)),
# Differencing = c(rep(d, times = model.length)))
# save prediction plot for preliminary review
arima.plot <- ggplot(arima.out, aes(x=Date, y = Hospitalizations_Total))+
geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), fill = "red", alpha = 0.5, size = 0.1)+
geom_line(color = "black", size = 1)+
labs(y = 'Daily Cases (7-Day Moving Average)', x = 'Date',
title = mydata$Level_Name[1]) +
scale_x_date(limits = c(mindate, maxdate + prediction.period),
date_labels = '%b-%d', date_breaks = '1 week') +
ggpubr::theme_pubr() +
theme(axis.text.x = element_text(angle = -90))
ggsave(plot = arima.plot, paste0('statistical-output/time-series/plots/',
mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'))
} else {
# insufficient data catch: return NA values for predictions
arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
# Cases_Raw = c(mydata$Hospitalizations_Total, rep(NA, times = prediction.period)),
Hospitalizations_Total = rep(NA, times = model.length),
CI_Lower = rep(NA, times = model.length),
CI_Upper = rep(NA, times = model.length))
# Order_AutoReg = c(rep(NA, times = model.length)),
# Order_Moving_Avg = c(rep(NA, times = model.length)),
# Differencing = c(rep(NA, times= model.length)))
}
#replace CI lower limit with 0 if negative
arima.out$CI_Lower <- ifelse(arima.out$CI_Lower>=0,arima.out$CI_Lower, 0)
return(arima.out)
}TSA - Hospitalization
# All data
TSA.hosp.arima.output.1 <- nlme::gapply(TSA.Daily,
FUN = covid.arima.forecast,
groups = TSA.Daily$Level_Name,
mindate = as.Date('2020-03-04'))## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
## Saving 7 x 5 in image
## Saving 7 x 5 in image## Warning: Removed 45 rows containing missing values (geom_path).
TSA.hosp.arima.df.1 <- data.table::rbindlist(TSA.hosp.arima.output.1, idcol=TRUE)
colnames(TSA.hosp.arima.df.1)[1] <- 'TSA'
TSA.hosp.arima.df.1 <- merge(TSA.hosp.arima.df.1, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))
# combine TSA
TSA.hosp.arima.df.1$TSA <- paste0(TSA.hosp.arima.df.1$TSA, ' - ', TSA.hosp.arima.df.1$TSA_Name)
TSA.hosp.arima.df.1$TSA_Name <- NULLState - hospitalization
texas.hosp.arima.df.1 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-03-04'))## Saving 7 x 5 in image
## Saving 7 x 5 in image
# texas.hosp.arima.df.2 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-06-03'))Hospitalization Write to CSV
tsa.hosp.arima<-write.csv(TSA.hosp.arima.df.1, "statistical-output/time-series/tsa_hosp_timeseries.csv", row.names=FALSE)
texas.hosp.arima<-write.csv(texas.hosp.arima.df.1, "statistical-output/time-series/state_hosp_timeseries.csv", row.names=FALSE)Stacking
TSA_hosp_out = TSA.hosp.arima.df.1
colnames(TSA_hosp_out)[1] = 'Level'
texas.hosp.arima.df.1$Level = 'Texas'
TSA_hosp_out$Level_Type = 'TSA'
texas.hosp.arima.df.1$Level_Type = 'State'
combined.hosp.arima.df.1 <- rbind(TSA_hosp_out, texas.hosp.arima.df.1)
write.csv(combined.hosp.arima.df.1, 'statistical-output/time-series/stacked_hosp_timeseries.csv',
row.names = FALSE)
write.csv(combined.hosp.arima.df.1, 'tableau/stacked_hosp_timeseries.csv',
row.names = FALSE)Hospitalization Forecast Plots
#Create plot function, forecast.data is the dataframe containing true data
# and the forecast data. days.forecast is the number of days that are forecast
# in the dataframe (here we are forecasting 10 days, so it should be 10 unless we
# change this at some point down the road)
forecast.plot<-function(forecast.data, plot.title, y.label){
mindate<-as.POSIXct("2020-03-04")
maxdate<-as.POSIXct(max(forecast.data$Date))
ts.plot<-ggplot(aes(x=as.POSIXct(Date), y=Hospitalizations_Total), data=forecast.data)+
geom_ribbon(aes(ymin=CI_Lower, ymax=CI_Upper), fill="grey50", size=0.1)+
geom_line(color="blue", size=1)+
geom_line(aes(x=as.POSIXct(Date), y=CI_Lower), color="grey", size=0.1)+
geom_line(aes(x=as.POSIXct(Date),y=CI_Upper), color="grey", size=0.1)+
scale_x_datetime(limits = c(mindate, maxdate))+
xlab("Date")+
ylab(y.label)+
#Can use the following title if we are running using nlme for data frame
#ggtitle(paste(forecast.data$.id,": TS Daily Cases + 10 Day Forcast",sep=""))
ggtitle(plot.title)
ts.plot
}######### View Output Table & Graph for TSA Q ##########
library(kableExtra)
#subset TSA Q arima data
TSA_Q.arima.df<-TSA.hosp.arima.df.1%>%subset(TSA.hosp.arima.df.1$TSA == "Q - Houston")
#apply the plot function to TSA Q arima data frame
TSA_Q.arima.df%>%kable(caption="Greater Houston ARIMA - Daily Hosp")%>%kable_styling(full_width=FALSE)| TSA | Date | Hospitalizations_Total | CI_Lower | CI_Upper |
|---|---|---|---|---|
| Q - Houston | 2020-03-04 | NA | NA | NA |
| Q - Houston | 2020-03-05 | NA | NA | NA |
| Q - Houston | 2020-03-06 | NA | NA | NA |
| Q - Houston | 2020-03-07 | NA | NA | NA |
| Q - Houston | 2020-03-08 | NA | NA | NA |
| Q - Houston | 2020-03-09 | NA | NA | NA |
| Q - Houston | 2020-03-10 | NA | NA | NA |
| Q - Houston | 2020-03-11 | NA | NA | NA |
| Q - Houston | 2020-03-12 | NA | NA | NA |
| Q - Houston | 2020-03-13 | NA | NA | NA |
| Q - Houston | 2020-03-14 | NA | NA | NA |
| Q - Houston | 2020-03-15 | NA | NA | NA |
| Q - Houston | 2020-03-16 | NA | NA | NA |
| Q - Houston | 2020-03-17 | NA | NA | NA |
| Q - Houston | 2020-03-18 | NA | NA | NA |
| Q - Houston | 2020-03-19 | NA | NA | NA |
| Q - Houston | 2020-03-20 | NA | NA | NA |
| Q - Houston | 2020-03-21 | NA | NA | NA |
| Q - Houston | 2020-03-22 | NA | NA | NA |
| Q - Houston | 2020-03-23 | NA | NA | NA |
| Q - Houston | 2020-03-24 | NA | NA | NA |
| Q - Houston | 2020-03-25 | NA | NA | NA |
| Q - Houston | 2020-03-26 | NA | NA | NA |
| Q - Houston | 2020-03-27 | NA | NA | NA |
| Q - Houston | 2020-03-28 | NA | NA | NA |
| Q - Houston | 2020-03-29 | NA | NA | NA |
| Q - Houston | 2020-03-30 | NA | NA | NA |
| Q - Houston | 2020-03-31 | NA | NA | NA |
| Q - Houston | 2020-04-01 | NA | NA | NA |
| Q - Houston | 2020-04-02 | NA | NA | NA |
| Q - Houston | 2020-04-03 | NA | NA | NA |
| Q - Houston | 2020-04-04 | NA | NA | NA |
| Q - Houston | 2020-04-05 | NA | NA | NA |
| Q - Houston | 2020-04-06 | NA | NA | NA |
| Q - Houston | 2020-04-07 | NA | NA | NA |
| Q - Houston | 2020-04-08 | NA | NA | NA |
| Q - Houston | 2020-04-09 | NA | NA | NA |
| Q - Houston | 2020-04-10 | NA | NA | NA |
| Q - Houston | 2020-04-11 | NA | NA | NA |
| Q - Houston | 2020-04-12 | NA | NA | NA |
| Q - Houston | 2020-04-13 | NA | NA | NA |
| Q - Houston | 2020-04-14 | NA | NA | NA |
| Q - Houston | 2020-04-15 | NA | NA | NA |
| Q - Houston | 2020-04-16 | NA | NA | NA |
| Q - Houston | 2020-04-17 | NA | NA | NA |
| Q - Houston | 2020-04-18 | 42.71429 | NA | NA |
| Q - Houston | 2020-04-19 | 50.85714 | NA | NA |
| Q - Houston | 2020-04-20 | 49.14286 | NA | NA |
| Q - Houston | 2020-04-21 | 46.00000 | NA | NA |
| Q - Houston | 2020-04-22 | 51.85714 | NA | NA |
| Q - Houston | 2020-04-23 | 57.42857 | NA | NA |
| Q - Houston | 2020-04-24 | 62.85714 | NA | NA |
| Q - Houston | 2020-04-25 | 68.28571 | NA | NA |
| Q - Houston | 2020-04-26 | 73.28571 | NA | NA |
| Q - Houston | 2020-04-27 | 86.85714 | NA | NA |
| Q - Houston | 2020-04-28 | 101.00000 | NA | NA |
| Q - Houston | 2020-04-29 | 105.42857 | NA | NA |
| Q - Houston | 2020-04-30 | 109.85714 | NA | NA |
| Q - Houston | 2020-05-01 | 114.00000 | NA | NA |
| Q - Houston | 2020-05-02 | 118.00000 | NA | NA |
| Q - Houston | 2020-05-03 | 121.85714 | NA | NA |
| Q - Houston | 2020-05-04 | 125.71429 | NA | NA |
| Q - Houston | 2020-05-05 | 129.57143 | NA | NA |
| Q - Houston | 2020-05-06 | 129.57143 | NA | NA |
| Q - Houston | 2020-05-07 | 133.42857 | NA | NA |
| Q - Houston | 2020-05-08 | 137.14286 | NA | NA |
| Q - Houston | 2020-05-09 | 140.57143 | NA | NA |
| Q - Houston | 2020-05-10 | 143.42857 | NA | NA |
| Q - Houston | 2020-05-11 | 146.85714 | NA | NA |
| Q - Houston | 2020-05-12 | 143.42857 | NA | NA |
| Q - Houston | 2020-05-13 | 136.71429 | NA | NA |
| Q - Houston | 2020-05-14 | 139.14286 | NA | NA |
| Q - Houston | 2020-05-15 | 141.42857 | NA | NA |
| Q - Houston | 2020-05-16 | 144.14286 | NA | NA |
| Q - Houston | 2020-05-17 | 147.14286 | NA | NA |
| Q - Houston | 2020-05-18 | 149.57143 | NA | NA |
| Q - Houston | 2020-05-19 | 158.71429 | NA | NA |
| Q - Houston | 2020-05-20 | 159.85714 | NA | NA |
| Q - Houston | 2020-05-21 | 162.85714 | NA | NA |
| Q - Houston | 2020-05-22 | 166.14286 | NA | NA |
| Q - Houston | 2020-05-23 | 169.14286 | NA | NA |
| Q - Houston | 2020-05-24 | 172.14286 | NA | NA |
| Q - Houston | 2020-05-25 | 175.14286 | NA | NA |
| Q - Houston | 2020-05-26 | 173.28571 | NA | NA |
| Q - Houston | 2020-05-27 | 186.85714 | NA | NA |
| Q - Houston | 2020-05-28 | 170.14286 | NA | NA |
| Q - Houston | 2020-05-29 | 170.85714 | NA | NA |
| Q - Houston | 2020-05-30 | 173.00000 | NA | NA |
| Q - Houston | 2020-05-31 | 175.14286 | NA | NA |
| Q - Houston | 2020-06-01 | 177.14286 | NA | NA |
| Q - Houston | 2020-06-02 | 183.85714 | NA | NA |
| Q - Houston | 2020-06-03 | 172.14286 | NA | NA |
| Q - Houston | 2020-06-04 | 193.28571 | NA | NA |
| Q - Houston | 2020-06-05 | 196.71429 | NA | NA |
| Q - Houston | 2020-06-06 | 198.57143 | NA | NA |
| Q - Houston | 2020-06-07 | 200.57143 | NA | NA |
| Q - Houston | 2020-06-08 | 202.42857 | NA | NA |
| Q - Houston | 2020-06-09 | 204.28571 | NA | NA |
| Q - Houston | 2020-06-10 | 224.28571 | NA | NA |
| Q - Houston | 2020-06-11 | 226.42857 | NA | NA |
| Q - Houston | 2020-06-12 | 228.71429 | NA | NA |
| Q - Houston | 2020-06-13 | 231.28571 | NA | NA |
| Q - Houston | 2020-06-14 | 233.85714 | NA | NA |
| Q - Houston | 2020-06-15 | 237.14286 | NA | NA |
| Q - Houston | 2020-06-16 | 240.71429 | NA | NA |
| Q - Houston | 2020-06-17 | 244.42857 | NA | NA |
| Q - Houston | 2020-06-18 | 248.71429 | NA | NA |
| Q - Houston | 2020-06-19 | 252.71429 | NA | NA |
| Q - Houston | 2020-06-20 | 257.14286 | NA | NA |
| Q - Houston | 2020-06-21 | 261.85714 | NA | NA |
| Q - Houston | 2020-06-22 | 266.42857 | NA | NA |
| Q - Houston | 2020-06-23 | 271.42857 | NA | NA |
| Q - Houston | 2020-06-24 | 276.57143 | NA | NA |
| Q - Houston | 2020-06-25 | 281.85714 | NA | NA |
| Q - Houston | 2020-06-26 | 288.00000 | NA | NA |
| Q - Houston | 2020-06-27 | 294.14286 | NA | NA |
| Q - Houston | 2020-06-28 | 300.42857 | NA | NA |
| Q - Houston | 2020-06-29 | 307.00000 | NA | NA |
| Q - Houston | 2020-06-30 | 313.42857 | NA | NA |
| Q - Houston | 2020-07-01 | 320.14286 | NA | NA |
| Q - Houston | 2020-07-02 | 327.00000 | NA | NA |
| Q - Houston | 2020-07-03 | 333.85714 | NA | NA |
| Q - Houston | 2020-07-04 | 341.14286 | NA | NA |
| Q - Houston | 2020-07-05 | 348.57143 | NA | NA |
| Q - Houston | 2020-07-06 | 356.00000 | NA | NA |
| Q - Houston | 2020-07-07 | 363.85714 | NA | NA |
| Q - Houston | 2020-07-08 | 371.57143 | NA | NA |
| Q - Houston | 2020-07-09 | 379.42857 | NA | NA |
| Q - Houston | 2020-07-10 | 387.57143 | NA | NA |
| Q - Houston | 2020-07-11 | 395.57143 | NA | NA |
| Q - Houston | 2020-07-12 | 403.71429 | NA | NA |
| Q - Houston | 2020-07-13 | 409.71429 | NA | NA |
| Q - Houston | 2020-07-14 | 418.14286 | NA | NA |
| Q - Houston | 2020-07-15 | 427.00000 | NA | NA |
| Q - Houston | 2020-07-16 | 435.14286 | NA | NA |
| Q - Houston | 2020-07-17 | 443.00000 | NA | NA |
| Q - Houston | 2020-07-18 | 450.42857 | NA | NA |
| Q - Houston | 2020-07-19 | 446.71429 | NA | NA |
| Q - Houston | 2020-07-20 | 456.28571 | NA | NA |
| Q - Houston | 2020-07-21 | 463.28571 | NA | NA |
| Q - Houston | 2020-07-22 | 469.85714 | NA | NA |
| Q - Houston | 2020-07-23 | 456.28571 | NA | NA |
| Q - Houston | 2020-07-24 | 463.38919 | 458.0047 | 468.7737 |
| Q - Houston | 2020-07-25 | 466.04615 | 458.3897 | 473.7026 |
| Q - Houston | 2020-07-26 | 472.93882 | 463.0025 | 482.8752 |
| Q - Houston | 2020-07-27 | 475.79659 | 463.9605 | 487.6327 |
| Q - Houston | 2020-07-28 | 482.49797 | 468.6002 | 496.3957 |
| Q - Houston | 2020-07-29 | 485.53796 | 469.7845 | 501.2914 |
| Q - Houston | 2020-07-30 | 492.06576 | 474.2892 | 509.8423 |
| Q - Houston | 2020-07-31 | 495.27110 | 475.6088 | 514.9334 |
| Q - Houston | 2020-08-01 | 501.64139 | 479.9374 | 523.3454 |
| Q - Houston | 2020-08-02 | 504.99677 | 481.3503 | 528.6432 |
#nlme::gapply(TSA.hosp.arima.df, FUN = forecast.plot, groups=TSA.daily$Level_Name)
forecast.plot(TSA_Q.arima.df, "TSA Q - Greater Houston Daily Hosp", "Daily Hosp")## Warning: Removed 45 rows containing missing values (geom_path).## Warning: Removed 142 rows containing missing values (geom_path).
## Warning: Removed 142 rows containing missing values (geom_path).
######### View Output Table & Graph for Texas ##########
texas.hosp.arima.df.1%>%kable(caption="Texas Arima - Daily Hosp")%>%kable_styling(full_width=FALSE)| Date | Hospitalizations_Total | CI_Lower | CI_Upper | Level | Level_Type |
|---|---|---|---|---|---|
| 2020-03-04 | 0.00000 | NA | NA | Texas | State |
| 2020-03-05 | 0.00000 | NA | NA | Texas | State |
| 2020-03-06 | 0.00000 | NA | NA | Texas | State |
| 2020-03-07 | 0.00000 | NA | NA | Texas | State |
| 2020-03-08 | 0.00000 | NA | NA | Texas | State |
| 2020-03-09 | 0.00000 | NA | NA | Texas | State |
| 2020-03-10 | 0.00000 | NA | NA | Texas | State |
| 2020-03-11 | 0.00000 | NA | NA | Texas | State |
| 2020-03-12 | 0.00000 | NA | NA | Texas | State |
| 2020-03-13 | 0.00000 | NA | NA | Texas | State |
| 2020-03-14 | 0.00000 | NA | NA | Texas | State |
| 2020-03-15 | 0.00000 | NA | NA | Texas | State |
| 2020-03-16 | 0.00000 | NA | NA | Texas | State |
| 2020-03-17 | 0.00000 | NA | NA | Texas | State |
| 2020-03-18 | 0.00000 | NA | NA | Texas | State |
| 2020-03-19 | 0.00000 | NA | NA | Texas | State |
| 2020-03-20 | 0.00000 | NA | NA | Texas | State |
| 2020-03-21 | 0.00000 | NA | NA | Texas | State |
| 2020-03-22 | 0.00000 | NA | NA | Texas | State |
| 2020-03-23 | 0.00000 | NA | NA | Texas | State |
| 2020-03-24 | 0.00000 | NA | NA | Texas | State |
| 2020-03-25 | 0.00000 | NA | NA | Texas | State |
| 2020-03-26 | 0.00000 | NA | NA | Texas | State |
| 2020-03-27 | 0.00000 | NA | NA | Texas | State |
| 2020-03-28 | 0.00000 | NA | NA | Texas | State |
| 2020-03-29 | 0.00000 | NA | NA | Texas | State |
| 2020-03-30 | 0.00000 | NA | NA | Texas | State |
| 2020-03-31 | 0.00000 | NA | NA | Texas | State |
| 2020-04-01 | 0.00000 | NA | NA | Texas | State |
| 2020-04-02 | 0.00000 | NA | NA | Texas | State |
| 2020-04-03 | 0.00000 | NA | NA | Texas | State |
| 2020-04-04 | 0.00000 | NA | NA | Texas | State |
| 2020-04-05 | 0.00000 | NA | NA | Texas | State |
| 2020-04-06 | 0.00000 | NA | NA | Texas | State |
| 2020-04-07 | 0.00000 | NA | NA | Texas | State |
| 2020-04-08 | 0.00000 | NA | NA | Texas | State |
| 2020-04-09 | 0.00000 | NA | NA | Texas | State |
| 2020-04-10 | 0.00000 | NA | NA | Texas | State |
| 2020-04-11 | 0.00000 | NA | NA | Texas | State |
| 2020-04-12 | 32.71429 | NA | NA | Texas | State |
| 2020-04-13 | 95.42857 | NA | NA | Texas | State |
| 2020-04-14 | 173.57143 | NA | NA | Texas | State |
| 2020-04-15 | 256.57143 | NA | NA | Texas | State |
| 2020-04-16 | 338.00000 | NA | NA | Texas | State |
| 2020-04-17 | 434.14286 | NA | NA | Texas | State |
| 2020-04-18 | 552.28571 | NA | NA | Texas | State |
| 2020-04-19 | 620.14286 | NA | NA | Texas | State |
| 2020-04-20 | 663.00000 | NA | NA | Texas | State |
| 2020-04-21 | 715.71429 | NA | NA | Texas | State |
| 2020-04-22 | 792.28571 | NA | NA | Texas | State |
| 2020-04-23 | 872.14286 | NA | NA | Texas | State |
| 2020-04-24 | 932.00000 | NA | NA | Texas | State |
| 2020-04-25 | 972.42857 | NA | NA | Texas | State |
| 2020-04-26 | 1025.57143 | NA | NA | Texas | State |
| 2020-04-27 | 1087.42857 | NA | NA | Texas | State |
| 2020-04-28 | 1117.85714 | NA | NA | Texas | State |
| 2020-04-29 | 1140.85714 | NA | NA | Texas | State |
| 2020-04-30 | 1143.85714 | NA | NA | Texas | State |
| 2020-05-01 | 1171.14286 | NA | NA | Texas | State |
| 2020-05-02 | 1168.42857 | NA | NA | Texas | State |
| 2020-05-03 | 1175.42857 | NA | NA | Texas | State |
| 2020-05-04 | 1167.14286 | NA | NA | Texas | State |
| 2020-05-05 | 1185.57143 | NA | NA | Texas | State |
| 2020-05-06 | 1191.85714 | NA | NA | Texas | State |
| 2020-05-07 | 1183.14286 | NA | NA | Texas | State |
| 2020-05-08 | 1189.42857 | NA | NA | Texas | State |
| 2020-05-09 | 1265.42857 | NA | NA | Texas | State |
| 2020-05-10 | 1295.14286 | NA | NA | Texas | State |
| 2020-05-11 | 1315.42857 | NA | NA | Texas | State |
| 2020-05-12 | 1321.85714 | NA | NA | Texas | State |
| 2020-05-13 | 1334.42857 | NA | NA | Texas | State |
| 2020-05-14 | 1329.71429 | NA | NA | Texas | State |
| 2020-05-15 | 1361.57143 | NA | NA | Texas | State |
| 2020-05-16 | 1307.57143 | NA | NA | Texas | State |
| 2020-05-17 | 1316.00000 | NA | NA | Texas | State |
| 2020-05-18 | 1310.71429 | NA | NA | Texas | State |
| 2020-05-19 | 1317.85714 | NA | NA | Texas | State |
| 2020-05-20 | 1346.85714 | NA | NA | Texas | State |
| 2020-05-21 | 1411.71429 | NA | NA | Texas | State |
| 2020-05-22 | 1396.85714 | NA | NA | Texas | State |
| 2020-05-23 | 1429.28571 | NA | NA | Texas | State |
| 2020-05-24 | 1466.42857 | NA | NA | Texas | State |
| 2020-05-25 | 1475.85714 | NA | NA | Texas | State |
| 2020-05-26 | 1461.57143 | NA | NA | Texas | State |
| 2020-05-27 | 1416.14286 | NA | NA | Texas | State |
| 2020-05-28 | 1414.28571 | NA | NA | Texas | State |
| 2020-05-29 | 1425.57143 | NA | NA | Texas | State |
| 2020-05-30 | 1454.00000 | NA | NA | Texas | State |
| 2020-05-31 | 1438.00000 | NA | NA | Texas | State |
| 2020-06-01 | 1440.42857 | NA | NA | Texas | State |
| 2020-06-02 | 1498.00000 | NA | NA | Texas | State |
| 2020-06-03 | 1488.71429 | NA | NA | Texas | State |
| 2020-06-04 | 1475.00000 | NA | NA | Texas | State |
| 2020-06-05 | 1472.85714 | NA | NA | Texas | State |
| 2020-06-06 | 1460.42857 | NA | NA | Texas | State |
| 2020-06-07 | 1433.00000 | NA | NA | Texas | State |
| 2020-06-08 | 1452.57143 | NA | NA | Texas | State |
| 2020-06-09 | 1477.00000 | NA | NA | Texas | State |
| 2020-06-10 | 1492.00000 | NA | NA | Texas | State |
| 2020-06-11 | 1510.85714 | NA | NA | Texas | State |
| 2020-06-12 | 1547.00000 | NA | NA | Texas | State |
| 2020-06-13 | 1564.57143 | NA | NA | Texas | State |
| 2020-06-14 | 1608.14286 | NA | NA | Texas | State |
| 2020-06-15 | 1661.42857 | NA | NA | Texas | State |
| 2020-06-16 | 1658.28571 | NA | NA | Texas | State |
| 2020-06-17 | 1763.42857 | NA | NA | Texas | State |
| 2020-06-18 | 1883.28571 | NA | NA | Texas | State |
| 2020-06-19 | 1948.57143 | NA | NA | Texas | State |
| 2020-06-20 | 2048.71429 | NA | NA | Texas | State |
| 2020-06-21 | 2149.85714 | NA | NA | Texas | State |
| 2020-06-22 | 2227.85714 | NA | NA | Texas | State |
| 2020-06-23 | 2303.71429 | NA | NA | Texas | State |
| 2020-06-24 | 2402.42857 | NA | NA | Texas | State |
| 2020-06-25 | 2437.14286 | NA | NA | Texas | State |
| 2020-06-26 | 2570.00000 | NA | NA | Texas | State |
| 2020-06-27 | 2642.85714 | NA | NA | Texas | State |
| 2020-06-28 | 2760.42857 | NA | NA | Texas | State |
| 2020-06-29 | 2810.85714 | NA | NA | Texas | State |
| 2020-06-30 | 2982.85714 | NA | NA | Texas | State |
| 2020-07-01 | 3088.57143 | NA | NA | Texas | State |
| 2020-07-02 | 3227.28571 | NA | NA | Texas | State |
| 2020-07-03 | 3370.00000 | NA | NA | Texas | State |
| 2020-07-04 | 3557.71429 | NA | NA | Texas | State |
| 2020-07-05 | 3649.85714 | NA | NA | Texas | State |
| 2020-07-06 | 3862.28571 | NA | NA | Texas | State |
| 2020-07-07 | 3994.00000 | NA | NA | Texas | State |
| 2020-07-08 | 4112.71429 | NA | NA | Texas | State |
| 2020-07-09 | 4308.00000 | NA | NA | Texas | State |
| 2020-07-10 | 4473.57143 | NA | NA | Texas | State |
| 2020-07-11 | 4510.42857 | NA | NA | Texas | State |
| 2020-07-12 | 4736.00000 | NA | NA | Texas | State |
| 2020-07-13 | 5023.57143 | NA | NA | Texas | State |
| 2020-07-14 | 5260.00000 | NA | NA | Texas | State |
| 2020-07-15 | 5328.14286 | NA | NA | Texas | State |
| 2020-07-16 | 5409.14286 | NA | NA | Texas | State |
| 2020-07-17 | 5348.14286 | NA | NA | Texas | State |
| 2020-07-18 | 5512.71429 | NA | NA | Texas | State |
| 2020-07-19 | 5565.28571 | NA | NA | Texas | State |
| 2020-07-20 | 5547.57143 | NA | NA | Texas | State |
| 2020-07-21 | 5540.57143 | NA | NA | Texas | State |
| 2020-07-22 | 5731.85714 | NA | NA | Texas | State |
| 2020-07-23 | 5740.14286 | NA | NA | Texas | State |
| 2020-07-24 | 5801.32587 | 5755.294 | 5847.358 | Texas | State |
| 2020-07-25 | 5862.50888 | 5786.979 | 5938.039 | Texas | State |
| 2020-07-26 | 5923.69189 | 5818.140 | 6029.244 | Texas | State |
| 2020-07-27 | 5984.87490 | 5847.718 | 6122.031 | Texas | State |
| 2020-07-28 | 6046.05791 | 5875.469 | 6216.646 | Texas | State |
| 2020-07-29 | 6107.24092 | 5901.351 | 6313.130 | Texas | State |
| 2020-07-30 | 6168.42393 | 5925.390 | 6411.458 | Texas | State |
| 2020-07-31 | 6229.60694 | 5947.632 | 6511.582 | Texas | State |
| 2020-08-01 | 6290.78995 | 5968.135 | 6613.445 | Texas | State |
| 2020-08-02 | 6351.97296 | 5986.954 | 6716.992 | Texas | State |
forecast.plot(texas.hosp.arima.df.1, "Texas Daily Hospitalizations", "Daily Hosp")## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 142 rows containing missing values (geom_path).
## Warning: Removed 142 rows containing missing values (geom_path).
STANDARD STATISTICAL TESTS
CASE RATIOS
TSA
#TODO: update all varnames
cumcases.TSA <- TSA.Daily %>% dplyr::select(Date,Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.TSA$Date <- ymd(cumcases.TSA$Date)
latestdate <- max(cumcases.TSA$Date)
earliestdate <- latestdate - 14
middate <- latestdate - 7
week2 <- subset(cumcases.TSA, Date==latestdate, select=Cases_Cumulative) -
subset(cumcases.TSA, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.TSA, Date==middate, select=Cases_Cumulative) -
subset(cumcases.TSA, Date==earliestdate, select=Cases_Cumulative)
current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
tsa_current_ratio_dat <- data.frame(TSA=subset(cumcases.TSA, Date==latestdate, select=Level_Name)[,1],
current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))
# add TSA name
tsa_current_ratio_dat <- merge(tsa_current_ratio_dat, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = 'TSA')
# combine TSA name
tsa_current_ratio_dat$TSA <- paste0(tsa_current_ratio_dat$TSA, ' - ', tsa_current_ratio_dat$TSA_Name)
tsa_current_ratio_dat$TSA_Name <- NULLPHR
cumcases.PHR <- PHR.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.PHR$Date <- ymd(cumcases.PHR$Date)
latestdate <- max(cumcases.PHR$Date)
earliestdate <- latestdate - 14
middate <- latestdate - 7
week2 <- subset(cumcases.PHR, Date==latestdate, select=Cases_Cumulative) -
subset(cumcases.PHR, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.PHR, Date==middate, select=Cases_Cumulative) -
subset(cumcases.PHR, Date==earliestdate, select=Cases_Cumulative)
current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
phr_current_ratio_dat <- data.frame(Level_Name=subset(cumcases.PHR, Date==latestdate, select=Level_Name)[,1],
current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))
# add PHR name
phr_current_ratio_dat <- merge(phr_current_ratio_dat, unique(PHR.Daily[, c('PHR', 'Level_Name')]), by = 'Level_Name')
# combine PHR name
phr_current_ratio_dat$Level_Name <- paste0(phr_current_ratio_dat$PHR, ' - ', phr_current_ratio_dat$Level_Name)
phr_current_ratio_dat$PHR <- NULL
colnames(phr_current_ratio_dat)[1] <- 'PHR'County
# declare vals
cumcases.county <- county.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.county$Date <- ymd(cumcases.county$Date)
latestdate <- max(cumcases.county$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7
# calc cumulative case ratios from 2 weeks ago and last week
week2 <- subset(cumcases.county, Date==latestdate, select=Cases_Cumulative) -
subset(cumcases.county, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.county, Date==middate, select=Cases_Cumulative) -
subset(cumcases.county, Date==earliestdate, select=Cases_Cumulative)
current_ratio <- (week2/week1)[,1]
# add contingencies
current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
# output
# TODO: investigate select vs subset(...)[, 'County']
county_current_ratio_dat <-
data.frame(County=subset(cumcases.county, Date==latestdate, select=Level_Name)[,1],
current_ratio=current_ratio,
current_ratio_cat = cut(current_ratio,
breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))Metro
cumcases.metro <- metro.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.metro$Date <- ymd(cumcases.metro$Date)
latestdate <- max(cumcases.metro$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7
week2 <- subset(cumcases.metro, Date==latestdate, select=Cases_Cumulative) -
subset(cumcases.metro, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.metro, Date==middate, select=Cases_Cumulative) -
subset(cumcases.metro, Date==earliestdate, select=Cases_Cumulative)
current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA
metro_current_ratio_dat <-
data.frame(Metro_Area=subset(cumcases.metro, Date==latestdate, select=Level_Name)[,1],
current_ratio=current_ratio,
current_ratio_cat = cut(current_ratio,
breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))State
cumcases.state <- state.Daily %>% dplyr::select(Date, Cases_Cumulative)
cumcases.state$Date <- ymd(cumcases.state$Date)
latestdate <- max(cumcases.state$Date, na.rm = TRUE)
earliestdate <- latestdate - 14
middate <- latestdate-7
week2 <- subset(cumcases.state, Date == latestdate, select = Cases_Cumulative) -
subset(cumcases.state, Date == middate, select = Cases_Cumulative)
week1 <- subset(cumcases.state, Date == middate, select = Cases_Cumulative) -
subset(cumcases.state, Date == earliestdate, select = Cases_Cumulative)
current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA
state_current_ratio_dat <-
data.frame(current_ratio=current_ratio,
current_ratio_cat = cut(current_ratio,
breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))export
write.csv(state_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'statistical-output/standard-stats/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'statistical-output/standard-stats/case-ratios/metro_case_ratio.csv', row.names = F)
write.csv(phr_current_ratio_dat,'statistical-output/standard-stats/case-ratios/phr_case_ratio.csv', row.names = F)
write.csv(state_current_ratio_dat, 'tableau/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'tableau/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'tableau/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'tableau/case-ratios/metro_case_ratio.csv', row.names = F)
write.csv(phr_current_ratio_dat,'tableau/case-ratios/phr_case_ratio.csv', row.names = F)grouping
colnames(county_current_ratio_dat)[1] <- 'Level'
colnames(metro_current_ratio_dat)[1] <- 'Level'
colnames(tsa_current_ratio_dat)[1] <- 'Level'
colnames(phr_current_ratio_dat)[1] <- 'Level'
state_current_ratio_dat$Level <- 'Texas'
county_current_ratio_dat$Level_Type = 'County'
metro_current_ratio_dat$Level_Type = 'Metro'
tsa_current_ratio_dat$Level_Type = 'TSA'
phr_current_ratio_dat$Level_Type = 'PHR'
state_current_ratio_dat$Level_Type = 'State'
combined_current_ratio_dat <- rbind(county_current_ratio_dat, metro_current_ratio_dat,
tsa_current_ratio_dat, state_current_ratio_dat, phr_current_ratio_dat)
stacked_ratio_out = combined_current_ratio_dat %>%
mutate(County = ifelse(Level_Type == 'County', as.character(Level), '')) %>%
mutate(TSA = ifelse(Level_Type == 'TSA', as.character(Level), '')) %>%
mutate(PHR = ifelse(Level_Type == 'PHR', as.character(Level), '')) %>%
mutate(Metro = ifelse(Level_Type == 'Metro', as.character(Level), '')) %>%
mutate(State = ifelse(Level_Type == 'State', as.character(Level), ''))
write.csv(stacked_ratio_out, 'statistical-output/standard-stats/case-ratios/stacked_case_ratio.csv',
row.names = FALSE)
# write.csv(combined_current_ratio_dat, 'tableau/stacked_case_ratio.csv',
# row.names = FALSE)% CHANGE
create.case.test <- function(level, dat, region){
# creates the % difference in cases and tests and smooth line with CIs
# level: either "TSA", "county", or "metro". Note that "county" won't work for many counties unless have enough cases.
# dat: dataset (e.g. "county", "metro", "tsa")
# region: the region within the dataset (county, metro region, or tsa)
if(level == "TSA"){
dat <- subset(dat, TSA==region)
}
if(level == "PHR"){
dat <- subset(dat, PHR==region)
}
if(level == "county"){
dat <- subset(dat, County == region)
}
if(level == "metro"){
dat <- subset(dat, Metro_Area==region)
}
# restrict data to first test date (for % test increase)
first.test.date <- ymd("2020-04-22")
dat$Date <- ymd(dat$Date)
dat$Date2 <- as.numeric(ymd(dat$Date))
# get 14 day rolling avg for instances where initial cases or tests are 0
dat$cases_avg14 <- rollmean(dat$Cases_Daily, k=14, fill = NA, align = 'right', na.rm = TRUE)
dat$tests_avg14 <- rollmean(dat$Tests_Daily, k=14, fill = NA, align = 'right', na.rm = TRUE)
# restrict new df to first.test.date forward
slopedata.tests <- data.frame(subset(dat, select=c(Date, Date2, Cases_Daily,
Tests_Daily, cases_avg14, tests_avg14))) %>%
subset(ymd(Date) >= ymd(first.test.date))
tests_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=Tests_Daily))
cases_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=Cases_Daily))
# numeric date for gam
tmpdata <- data.frame(Date2=slopedata.tests$Date2)
if (cases_start != 0 & tests_start != 0) {
# Calc splines if cases and tests will produce non Inf values
slopedata.tests$newtestsY <- 100*(as.vector(slopedata.tests$Tests_Daily / tests_start) - 1)
slopedata.tests$newcasesY <- 100*(as.vector(slopedata.tests$Cases_Daily / cases_start) - 1)
# fit and predict w/ gam for spline vals
cases.gam <- gam(newcasesY~s(Date2,bs="cs", k=5), data=slopedata.tests)
tests.gam <- gam(newtestsY~s(Date2,bs="cs",k=5), data=slopedata.tests)
cases.fit <- predict(cases.gam, tmpdata, se.fit=TRUE)
tests.fit <- predict(tests.gam, tmpdata, se.fit=TRUE)
# Use 95% CI for splines
tmp.df <- data.frame(Date = slopedata.tests$Date,
Data_Type = 'Raw',
cases_percentdiff = slopedata.tests$newcasesY,
tests_percentdiff = slopedata.tests$newtestsY,
cases_percentdiff_spline = cases.fit$fit,
cases_percentdiff_spline_lower = cases.fit$fit-1.96*cases.fit$se,
cases_percentdiff_spline_upper = cases.fit$fit+1.96*cases.fit$se,
tests_percentdiff_spline = tests.fit$fit,
tests_percentdiff_spline_lower = tests.fit$fit-1.96*tests.fit$se,
tests_percentdiff_spline_upper = tests.fit$fit+1.96*tests.fit$se)
} else {
# only calc % increase for regions w/ sparse cases
cases_avg_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=cases_avg14))
tests_avg_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=tests_avg14))
slopedata.tests$newcasesY <- 100*(as.vector(slopedata.tests$cases_avg14 / cases_avg_start) - 1)
slopedata.tests$newtestsY <- 100*(as.vector(slopedata.tests$tests_avg14 / tests_avg_start) - 1)
tmp.df <- data.frame(Date = slopedata.tests$Date,
Data_Type = '14-Day Average',
cases_percentdiff = slopedata.tests$newcasesY,
tests_percentdiff = slopedata.tests$newtestsY,
cases_percentdiff_spline = rep(NA, times = nrow(slopedata.tests)),
cases_percentdiff_spline_lower = rep(NA, times = nrow(slopedata.tests)),
cases_percentdiff_spline_upper = rep(NA, times = nrow(slopedata.tests)),
tests_percentdiff_spline = rep(NA, times = nrow(slopedata.tests)),
tests_percentdiff_spline_lower = rep(NA, times = nrow(slopedata.tests)),
tests_percentdiff_spline_upper = rep(NA, times = nrow(slopedata.tests)))
}
return(tmp.df)
}State
state <- readxl::read_xlsx('combined-datasets/state.xlsx', sheet=1)
state.case.test.df <- create.case.test(level="State", state, NA)
write.csv(state.case.test.df, 'statistical-output/standard-stats/pct-change/state_pct_change.csv',
row.names = FALSE)TSA
tsa <- read.csv('combined-datasets/tsa.csv')
# Create tsa-level data (can optimize this code in the future, but it runs pretty quickly already)
# TODO: convert to gapply -> rbindlist
tsalist <- unique(tsa$TSA)
tmp <- create.case.test(level="TSA", tsa,tsalist[1])
tsa.case.test.df <- data.frame(TSA=rep(tsalist[1], nrow(tmp)),
TSA_Name=rep(unique(tsa$TSA_Name[1])),
tmp)
for(i in c(2:length(tsalist))){
tmp <- create.case.test(level="TSA", tsa,tsalist[i])
tsa.case.test.df <- rbind(tsa.case.test.df, data.frame(TSA=rep(tsalist[i], nrow(tmp)),
TSA_Name=rep(unique(tsa$TSA_Name[i])),
tmp))
}
## combine tsa name
tsa.case.test.df$TSA = paste0(tsa.case.test.df$TSA, ' - ', tsa.case.test.df$TSA_Name)
tsa.case.test.df$TSA_Name = NULL
write.csv(tsa.case.test.df, 'statistical-output/standard-stats/pct-change/tsa_pct_change.csv',
row.names = FALSE)PHR
phr <- read.csv('combined-datasets/phr.csv')
# Create phr-level data (can optimize this code in the future, but it runs pretty quickly already)
# TODO: convert to gapply -> rbindlist
phrlist <- unique(phr$PHR)
tmp <- create.case.test(level="PHR", phr,phrlist[1])
phr.case.test.df <- data.frame(PHR=rep(phrlist[1], nrow(tmp)),
PHR_Name=rep(unique(phr$PHR_Name[1])),
tmp)
for(i in c(2:length(phrlist))){
tmp <- create.case.test(level="PHR", phr,phrlist[i])
phr.case.test.df <- rbind(phr.case.test.df, data.frame(PHR=rep(phrlist[i], nrow(tmp)),
PHR_Name=rep(unique(phr$PHR_Name[i])),
tmp))
}
## combine phr name
phr.case.test.df$PHR = paste0(phr.case.test.df$PHR, ' - ', phr.case.test.df$PHR_Name)
phr.case.test.df$PHR_Name = NULL
write.csv(phr.case.test.df, 'statistical-output/standard-stats/pct-change/phr_pct_change.csv',
row.names = FALSE)Metro
Metro <- read.csv('combined-datasets/metro.csv')
metrolist <- unique(Metro$Metro_Area)
tmp <- create.case.test(level="metro", Metro,metrolist[1])
metro.case.test.df <- data.frame(Metro_Area=rep(metrolist[1], nrow(tmp)), tmp)
tmp <- create.case.test(level="metro", Metro,metrolist[2])
metro.case.test.df <- rbind(metro.case.test.df, data.frame(Metro_Area=rep(metrolist[2], nrow(tmp)), tmp))
write.csv(metro.case.test.df, 'statistical-output/standard-stats/pct-change/metro_pct_change.csv',
row.names = FALSE)County
county <- read.csv('combined-datasets/county.csv')
countylist <- unique(county$County)
tmp <- create.case.test(level="county", county,countylist[1])
county.case.test.df <- data.frame(County=rep(countylist[1], nrow(tmp)), tmp)
for(i in 2:length(countylist)) {
tmp <- create.case.test(level="county", county,countylist[i])
county.case.test.df <- rbind(county.case.test.df, data.frame(County=rep(countylist[i], nrow(tmp)), tmp))
}
# assess missing vals
missing_vals <- sum(is.na(county.case.test.df$cases_percentdiff)) + sum(is.na(county.case.test.df$tests_percentdiff))
recorded_vals <- sum(!is.na(county.case.test.df$cases_percentdiff)) + sum(!is.na(county.case.test.df$tests_percentdiff))
# 11.88% w/ 2020-05-01
# 12.71% w/ 2020-04-22
# 12% w/ 2020-05-15
missing_vals / (missing_vals + recorded_vals)## [1] 0.0982347write.csv(county.case.test.df, 'statistical-output/standard-stats/pct-change/county_pct_change.csv',
row.names = FALSE)grouping
colnames(county.case.test.df)[1] <- 'Level'
colnames(metro.case.test.df)[1] <- 'Level'
colnames(tsa.case.test.df)[1] <- 'Level'
colnames(phr.case.test.df)[1] <- 'Level'
state.case.test.df$Level <- 'Texas'
county.case.test.df$Level_Type = 'County'
metro.case.test.df$Level_Type = 'Metro'
tsa.case.test.df$Level_Type = 'TSA'
phr.case.test.df$Level_Type = 'PHR'
state.case.test.df$Level_Type = 'State'
combined.case.test.df <- rbind(county.case.test.df, metro.case.test.df,
tsa.case.test.df, state.case.test.df, phr.case.test.df)
write.csv(combined.case.test.df, 'statistical-output/standard-stats/pct-change/stacked_pct_change.csv',
row.names = FALSE)
write.csv(combined.case.test.df, 'tableau/stacked_pct_change.csv',
row.names = FALSE)STACK COMBINATIONS
combined_current_ratio_dat$Date = max(combined.case.test.df$Date)
colnames(combined.arima.df.1)[5:6] = c('TS_CI_Lower', 'TS_CI_Upper')
colnames(combined.rt.df)[4:5] = c('RT_CI_Lower', 'RT_CI_Upper')
stacked_all = Reduce(function(x, y) merge(x, y, by = c('Level_Type', 'Level', 'Date'), all=TRUE),
list(combined.case.test.df, combined_current_ratio_dat,
combined.arima.df.1, combined.rt.df, combined.hosp.arima.df.1))
write.csv(stacked_all, 'tableau/stacked_critical_trends.csv', row.names = FALSE)rolling average
Omitted fow now, Dr. Yaseen is processing this in Tableau
Trends
Requires substantial cleaning & consideration for visualization - omitting until dashboard version 2